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Chapter  1 
Introduction 


The  potential  of  sustained  hypersonic  flight  to  revolutionize  military  and  commercial  activ¬ 
ity  is  well  recognized,  and  is  reflected  in  recent  initiatives  such  as  the  National  Aerospace 
Initiative.  High-speed  vehicles  will  substantially  impact  military  strategy  by  providing  new 
defensive  options  such  as  a  rapid  on-demand  global  strike  capability  with  much  shorter  re¬ 
sponse  times  than  currently  possible.  Furthermore,  the  development  of  new  technologies 
based  on  air-breathing  propulsion  can  be  leveraged  to  considerably  reduce  the  cost  of  access- 
to-space,  the  benehts  of  which  are  both  military  as  well  as  commercial. 

However,  daunting  technical  challenges  remain  in  realizing  such  vehicles.  The  harsh 
environment  imposed  by  the  envelope  of  such  future  missions  is  manifested  in  the  severe 
anticipated  thermo-mechanical  loads  and  various  propulsion-related  requirements.  Although 
the  diversity  of  the  physical  phenomena  encountered  is  broad,  several  key  limiting  issues  have 
been  identihed  as  primary  challenges,  including  both  local  and  global  constraints  such  as, 
for  example,  cowl  lip  loading  and  airframe  balance.  A  scrutiny  of  the  problems  identihed 
reveals  the  pervasive  importance  of  several  basic  huid  dynamic  phenomena.  One  of  these, 
and  possibly  the  least  understood,  is  that  of  high-speed  transition. 

In  the  present  report,  the  main  results  stemming  from  research  supported  by  AFOSR 
grant  are  presented.  The  theoretical  and  computational  studies  of  stability,  transition  and 
how  control  have  been  carried  out  with  an  emphasis  on  the  multimode  decomposition  in¬ 
cluding  nonparallel  how  ehects.  The  multimode  decomposition  can  serve  to  analysis  of  DNS 
results  for  transitional  boundary  layers.  Application  of  the  multimode  decomposition  re¬ 
quired  re-evaluation  of  discrete  modes  branching  and  clarihcation  of  the  terminology  issues. 
Ehect  of  chemical  reactions  in  binary  mixtures  of  oxygen  and  nitrogen  has  been  studied 
within  the  inviscid  stability  analysis.  In  addition  to  the  studies  of  high-speed  boundary 
layers,  fundamental  aspects  were  addressed  to  evaluate  fundamentals  of  global  (bi-global) 
stability  concept. 


The  principal  investigator  is  thankful  to  Prof.  A.  Fedorov,  Prof.  X.  Zhong,  Dr.  J. 
Klentzman,  Dr.  X.  Wang,  Mr.  E.  Ulker  for  their  signihcant  input  into  this  three-year 
project. 
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The  results  have  been  published  in  Encyclopedia  of  Aerospace  Engineering,  1  journal 
paper  and  in  5  conference  publications. 

Encyclopedia  of  Aerospace  Engineering 

1.  A.  Tumin,  Flow  Instabilities  and  Transition,  in  Encyclopedia  of  Aerospace  Engineering, 
R.  Blocley  and  W.  Shyy  (eds),  John  Wiley  &  Sons  Ltd,  Chichester,  UK,  pp.  139-150 

Journal  Paper 

1.  A.  Tumin,  X.  Wang,  and  X.  Zhong,  Numerical  Simulation  and  Theoretical  Analysis  of 
Perturbations  in  Hypersonic  Boundary  Layers,  AIAA  J,  Vol.  49,  No.  3,  2011. 

Conference  publications 

1.  A.  Tumin,  X.  Wang,  and  X.  Zhong,  Numerical  Simulation  and  Theoretical  Analysis  on 
Hypersonic  Boundary-Layer  Receptivity  to  Wall  Blowing-Suction,  Orlando  EL,  AIAA 
Paper  2010-0534,  2010. 

2.  A.  Fedorov  and  A.  Tumin,  Branching  of  Discrete  Modes  in  High-Speed  Boundary 
Layers  and  Terminology  Issues,  Chicago  IL,  AIAA  Paper  2010-5003,  2010. 

3.  E.  Ulker,  J.  Klentzman,  and  A.  Tumin,  Stability  of  Boundary  Layers  in  Binary  Mixtures 
of  Oxygen  and  Nitrogen,  Orlando  EL,  AIAA  Paper  2011-0370,  2011. 

4.  A.  Tumin,  X.  Wang,  and  X.  Zhong,  Direct  Numerical  Simulation  and  Theoretical  Anal¬ 
ysis  of  Perturbations  in  Hypersonic  Boundary  Layers,  in  Seventh  lUTAM  Symposium 
on  Laminar- Turbulent  Transition,  Proceedings  of  the  Seventh  lUTAM  Symposium  on 
Laminar-Turbulent  Transition,  Stockholm,  Sweden,  2009,  Springer,  2010,  pp.  427-432. 

5.  A.  Tumin,  Toward  the  Foundation  of  a  Global  (Bi-Global)  Modes  Concept,  in  Global 
Elow  Instability  and  Control-IV,  Greta  Maris,  Hersinissos,  Crete,  Sept  28-Oct  2,  2009, 
V.  Theohlis,  T.  Colonius,  and  A.  Seifert  eds.,  ISBN-13:  978-84-692-6247-4 
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Chapter  2 


Numerical  simulation  and  theoretical 
analysis  on  hypersonic  boundary- layer 
receptivity  to  wall  blowing-suction 


This  section  is  based  on  the  results  published  in  colaboration  with  X.  Wang  and  X.  Zhong 
[TWZIO,  TWZll] 

2.1  Introduction 

The  progress  being  made  in  computational  fluid  dynamics  provides  an  opportunity  for  reli¬ 
able  simulations  of  such  complex  phenomena  as  laminar-turbulent  transition.  The  dynam¬ 
ics  of  flow  transition  depends  on  the  instability  of  small  perturbations  excited  by  external 
sources.  Computational  results  provide  complete  information  about  the  flow  field  that  would 
be  impossible  to  measure  in  real  experiments. 

Recently,  a  method  of  normal  mode  decomposition  was  developed  for  two-  and  three- 
dimensional  perturbations  in  compressible  and  incompressible  boundary  layers  [Tum03,  GT04, 
Tum07].  In  Ref.  [TWZ07],  the  method  was  applied  to  the  theoretical  analysis  of  a  pertur¬ 
bation  flow  field  in  the  vicinity  of  a  blowing-suction  actuator  obtained  from  direct  numerical 
simulation  (DNS).  The  results  demonstrated  very  good  agreement  between  the  amplitudes  of 
the  modes  Altered  out  from  the  DNS  data  and  those  solved  by  linear  theory  of  the  flow  recep¬ 
tivity  to  wall  blowing-suction.  However,  the  development  of  the  perturbations  downstream 
from  the  actuator  has  not  been  analyzed  yet. 

Perturbations  observed  in  experiments  and  computations  in  the  vicinity  of  an  actuator 
possess  a  nonmonotonic  character.  This  behavior  occurs  because  the  perturbation  intro¬ 
duced  by  the  actuator  is  composed  of  modes  of  the  discrete  (unstable  and  stable  modes)  and 
continuous  spectra,  and  one  cannot  distinguish  the  unstable  mode  clearly.  Are  the  observa¬ 
tions  still  compatible  with  the  linear  stability  theory  (LST)?  In  order  to  answer  this  question, 
we  must  decompose  the  perturbation  into  the  normal  modes  and  compare  their  amplitudes 
with  those  predicted  by  LST.  However,  the  LST  prediction  must  take  into  account  the  non- 
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parallel  boundary- layer  flow  effects  because  the  development  of  the  perturbation  takes  place 
on  a  length  scale  much  larger  than  the  boundary- layer  thickness. 

The  nonparallel  flow  effects  on  the  development  of  unstable  discrete  modes  on  a  length 
scale  that  is  much  larger  than  the  boundary- layer  thickness  have  been  studied  within  the 
scope  of  the  method  of  multiple  scales  [Bou72,  Gas74,  SN75,  SN77,  PN79,  Gap80,  NaySO, 
EH80,  TF82]  (MMS).  Another  method  that  allows  inclusion  of  the  nonparallel  flow  effects 
on  unstable  modes  is  based  on  the  parabolized  stability  equations  [HB87,  Ber91,  Her97] 
(PSE).  Fedorov  and  Khokhlov  [FKOl]  pointed  out  that  the  role  of  decaying  modes  can  be 
significant,  and  one  has  must  pay  attention  to  them  if  there  is  a  synchronism  with  the 
other  modes.  Because  the  analysis  of  decaying  modes  is  important,  one  must  include  the 
nonparallel  boundary- layer  flow  effects  using  the  MMS. 

In  the  present  work,  we  apply  the  multimode  decomposition  to  DNS  results  downstream 
from  the  blowing-suction  actuator  in  hypersonic  boundary  layers  past  a  flat  plate  and  a 
sharp  wedge  to  compare  the  amplitudes  of  the  modes  found  from  the  computations  with  the 
prediction  of  the  linear  stability  theory  including  nonparallel  flow  effects. 


2.2  Outline  of  the  multimode  decomposition 


The  method  of  multimode  decomposition  of  perturbations  having  a  prescribed  frequency 
is  based  on  the  biorthogonal  eigenfunction  system  for  linearized  Navier-Stokes  equations 
[Tum07].  For  the  clarity  of  further  discussion,  we  reproduce  the  main  definitions  necessary 
for  discussing  the  present  work. 

We  consider  a  compressible  two-dimensional  boundary  layer  in  Gartesian  coordinates, 
where  x  and  2:  are  the  downstream  and  spanwise  coordinates,  respectively,  and  coordinate  y 
corresponds  to  the  distance  from  the  wall.  We  write  the  linearized  Navier-Stokes  equations 
for  a  periodic-in-time  perturbation  (the  frequency  is  equal  to  zero  in  the  case  of  a  roughness- 
induced  perturbation),  ~  exp  {—iut),  in  matrix  form  as 


d  /  dA\  dA  „  ^  ^  dA  ^  dA  „  ^ 

~a~  (  “  HiA  -|-  H2^ - h  H3— - h  H4A, 

oy  \  oy  J  oy  ox  oz 

where  vector  A  has  16  components 

A  (x,  y,  z)  =(u,  du/dy,  v,  tt,  9, 06 /dy,  w,  dw/dy,  du/dx,  dv/dx, 

d6/dx,  dw/dx,  du/dz,  dv/dz,  06 /Oz,  Ow/OzY' . 


(2.1) 


(2.2) 


Lo,Li,Hi,H2,H3,  and  H4  are  16  x  16  matrices  (their  definitions  are  given  in  Ref.  [Tum08]); 
u,  u,  tc,  TT,  and  6  represent  three  velocity  components,  pressure,  and  temperature  perturba¬ 
tions,  respectively;  and  the  superscript  T  in  (2.2)  stands  for  transpose.  Matrix  H4  originates 
from  the  nonparallel  character  of  the  flow.  It  includes  terms  with  the  y-component  of  the 
mean  flow  velocity  and  derivatives  of  the  mean  flow  profiles  with  respect  to  the  coordinate 

X. 

In  the  quasi-parallel  flow  approximation,  the  solution  of  the  linearized  Navier-Stokes 
equations  can  be  expanded  into  normal  modes  of  the  discrete  and  continuous  spectra  {Aq,^,  Bq,^} 
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[Tum07],  where  Aa/3  and  are  eigenfunctions  of  the  direct  and  adjoint  problems.  Sub¬ 
scripts  a  and  f3  indicate  the  eigenfunctions  corresponding  to  the  streamwise,  a,  and  spanwise, 
/5,  wave  numbers,  respectively.  The  eigenfunction  system  { has  an  orthogonality 
relation  given  as 


(H2Aa/3,  Bq,//3)  =  j  (H2Aq,/3,  —  TAq-q,/,  (2.3) 

0 

where  T  is  a  normalization  constant,  is  a  Kronecker  delta  if  either  a  or  a'  belongs  to 
the  discrete  spectrum,  and  Aqq,/  is  a  Dirac  delta  function  if  both  a  and  a'  belong  to  the 
continuous  spectrum. 

The  orthogonality  relation  (2.3)  can  be  used  for  decomposition  of  the  DNS  data  at  x  =  xq 
into  normal  modes.  After  applying  the  Fourier  transform  in  ;2,  the  data  can  be  presented  in 
the  form  of  a  vector 


A/30  {xq,  y)  =  ^  {P)  (xo,  y) ,  (2.4) 

Q; 

where  ^  stands  for  summation  and  integration  over  the  discrete  and  continuous  spectra, 
respectively.  Applying  the  orthogonality  relation  (2.3),  one  can  hnd  the  coefficients  Ca  {P) 
as  follows: 


cn;3)  =  (2,5) 

In  the  case  of  the  transient  growth  phenomenon  studied  in  experiments  [Whi02,  WE03, 
FBTC04,  WRE05,  DW08]  and  DNS  [FC04,  RV07,  SAD+09],  the  flow  perturbation  is  given 
only  by  modes  of  the  continuous  spectrum.  The  multimode  decomposition  has  been  applied 
by  Denissen  and  White  [DW09]  to  the  analysis  of  the  DNS  data  of  Rizzetta  and  Visbal 
[RV07].  The  authors  demonstrated  that  the  method  can  be  utilized  even  when  only  partial 
information  in  vector  A^o  {xo,y)  is  available. 

In  a  weakly  nonparallel  flow,  one  can  employ  the  method  of  multiple  scales  (MMS)  by 
introducing  fast  {x)  and  slow  (A  =  ex,  e  1)  scales.  The  mean  flow  prohles  depend  on  y 
and  X  only,  whereas  the  perturbation  will  depend  on  both  length  scales.  In  the  case  of  a 
discrete  mode,  the  solution  of  the  linearized  Navier-Stokes  equation  is  presented  in  the  form 


(i,  A,  t;)  =  (A')  A™,,  (A,  y)  +  eA™^  (A,  y)  + 


(2.6) 


where  the  function  (A)  must  be  determined.  After  substitution  of  Eq.  (2.6)  into  Eq. 
(2.1),  we  arrive  in  order  0{e)  at  an  inhomogeneous  equation  for  The  solvability 

condition  of  this  equation  allows  for  the  Ending  of  D^y  (A)  (details  and  relevant  references 
can  be  found  in  Ref.  [Tum08]). 
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2.3  Direct  numerical  simulation  approach 

In  direct  numerical  simulation,  the  receptivity  of  hypersonic  boundary  layers  over  a  flat  plate 
and  a  sharp  wedge  to  wall  blowing-suction  are  considered  by  solving  the  two-dimensional 
compressible  Navier-Stokes  equations.  Wall  blowing-suction  is  introduced  by  a  slot  located 
near  the  leading  edge.  In  the  assumption  of  thermally  and  calorically  perfect  gas  flows,  the 
governing  equations  in  conservative  variables  are  given  as 

where  is  a  column  vector  containing  the  conservative  variables, 

U  =  {p,pu,f)v,eY .  (2.8) 

The  flux  vectors  in  (2.7)  are  divided  into  their  inviscid  and  viscous  components  due  to  the 
fact  that  the  two  components  are  discretized  with  different  schemes.  The  components,  Fn 
and  7^2*,  are  inviscid  flux  whereas  Fiv  and  F2v  are  viscous  flux  components. 


pUj 

puuj  -|-  pAij 
pVUj  +  pA2j 
Uj{e  +  p) 


(2.9) 


0 


'^XXj 

—Tyxj 

'T’xh-Xj'^k  F 


dT 

dxj  J 


(2.10) 


with  j,  k  G  (1, 2).  In  Cartesian  coordinates,  Xi  and  Ui  are  dehned  in  the  streamwise  direction 
{x  and  u)  whereas  X2  and  U2  are  dehned  in  the  wall-normal  direction  {y  and  v). 

Under  the  perfect  gas  assumption,  pressure  and  energy  are  given  by 


p  =  pRT, 


(2.11) 


e  =  pc^T  +  +  v^),  (2.12) 

where  c„  is  the  specihc  heat  at  constant  volume.  In  the  simulation,  the  viscosity  coefficient,  p, 
and  the  heat  conductivity,  K,  are  calculated  using  Sutherland’s  law  together  with  a  constant 
Prandtl  number,  Pr.  They  are  both  functions  of  temperature  only. 


3/2 


Tr  +  Ts 

T  +  T,  ’ 


p  =  P' 


T 

Tr 


(2.13) 


(2.14) 


K 


fJjCp 


where  =  1-7894  x  10“^  Ns/m^,  =  288.0  K,  T*  =  110.33  K,  and  Cp  is  the  specific  heat 
at  constant  pressure. 

The  scales  of  velocity  and  length  are  the  free-stream  velocity  U^o  and  the  Blasius  length 
scale,  L  =  {jJioox/ PooUooY^"^ ■  Since  the  length  scale  is  changing  along  the  fiat  plate,  it  is  not 
convenient  to  scale  the  frequency  in  terms  of  L/Uoo-  In  stead,  the  frequency  is  scaled  as 
follows: 

p  _  27r/ /ipo  _  ujpoo  ,  . 

o  U‘^  o  1 

roQ^  oo  KOODOO 


In  boundary  layer  analysis,  the  Reynolds  number  based  on  the  Blasius  length  scale  is  gener¬ 
ally  used, 

r=P^U^L  ^2,16) 

1^00 

In  the  present  work,  the  fifth-order  shock-fitting  finite  difference  method  of  Zhong  [Zho98] 
is  used  to  solve  the  governing  equations  in  a  domain  bounded  by  the  bow  shock  and  the  fiat 
plate  (or  wedge).  In  other  words,  the  bow  shock  is  treated  as  a  boundary  of  the  computational 
domain.  The  Rankine-Hugoniot  relations  across  the  shock  and  a  characteristic  compatibility 
relation  coming  from  the  downstream  flow  field  are  combined  to  solve  for  the  flow  variables 
behind  the  shock.  The  shock-fitting  method  makes  it  possible  for  the  Navier-Stokes  equations 
to  be  spatially  discretized  by  high-order  finite  difference  methods.  Specifically,  a  fifth-order 
upwind  scheme  is  applied  to  discretize  the  inviscid  flux  derivatives.  By  using  the  shock¬ 
fitting  method,  the  interaction  between  the  bow  shock  and  the  wall  blowing-suction  induced 
perturbations  is  solved  as  a  part  of  the  solutions  with  the  position  and  velocity  of  the  shock 
front  being  solved  as  dependent  variables.  Both  of  the  cases  correspond  to  an  adiabatic  wall 
boundary  condition. 

The  numerical  simulations  in  the  current  paper  have  been  well  validated.  For  the  flow 
over  the  fiat  plate,  three  sets  of  grid  structures  are  used  to  check  the  grid  independence  of 
numerical  simulation  results  near  the  leading  edge.  The  steady  flow  is  compared  with  the 
experimental  measurements  of  Maslov  et  al.[MSSA01].  At  three  different  locations,  a;  =  96 
mm,  121  mm,  and  138  mm  {R  =  1134.46, 1254.19, 1329.66),  the  distributions  of  the  dimen¬ 
sionless  streamwise  velocity  and  normalized  Mach  number  in  the  wall-normal  direction  are 
in  good  agreement  with  those  measured  in  experiments.  Two  sets  of  the  grid  structures 
are  used  to  check  grid  independence  of  unsteady  numerical  simulations.  The  comparison  of 
pressure  perturbation  amplitudes  calculated  using  the  two  grid  structures  shows  that  the 
grid  structure  used  for  the  simulation  is  sufficient.  Details  of  the  validations  are  presented 
in  Ref.  [WZ09].  In  Ref.  [TWZ07],  validations  of  the  numerical  simulation  results  on  the 
sharp  wedge  were  carried  out  for  both  steady  and  unsteady  flows.  Furthermore,  the  nu¬ 
merical  perturbation  field  just  downstream  of  the  blowing-suction  slot  was  decomposed  into 
boundary-layer  waves.  The  amplitudes  of  the  decomposed  waves  were  in  good  agreement 
with  those  obtained  from  a  theoretical  receptivity  model. 
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Figure  2.1:  Amplitude  distribution  of  the  blowing-suction  through  the  slot. 


2.4  Results 


2.4.1  Flat  plate 

The  free-stream  flow  conditions  that  we  consider  are:  Mach  number  =  5.92,  temperature 
Too  =  48.69  K,  pressure  p^o  =  742.76  Pa.  The  Prandtl  number  and  the  specihc  heat  ratio  are 
0.72  and  1.4,  respectively.  The  free  stream  flow  parameters  are  the  same  as  those  of  Maslov 
et  al.  [MSSAOl]. 

The  viscosity  coefficient  is  calculated  by  Sutherland’s  law.  The  dimensionless  blowing- 
suction  of  mass  flux  at  the  wall  is  expressed  as 


(pv)' =  eg  (1)  S{t),  e  =  0.405  X  10 

_  f  20.25/5  -  35.4375/^  +  15.1875/2,  (/<!); 

®  ^  1  -20.25(2-  i)“  + 35.4375  (2  15.1875(2  (/>!), 


(2.17) 


2(x  —  Xi) 

1[X)  =  - - -,  Xi<X<Xe, 

{Xe  -  Xi) 


where  Xi  =  33  mm  and  Xe  =  37  mm  are  the  coordinates  of  the  leading  and  the  trailing  edges 
of  the  slot,  respectively.  The  amplitude  distribution,  g  (/),  is  shown  in  Fig.  2.1. 

The  function  of  time  S{t)  in  (2.17)  is  dehned  as 


S{t) 


1,  mod(/,  20/is)  <  2/is; 
0,  mod(/,  20/is)  >  2ps. 


(2.18) 
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Figure  2.2:  Discrete  modes  and  the  continuous  spectrum. 


The  function  S{t)  can  be  expressed  as  a  Fourier  series. 

Analyses  of  the  mean  flow  velocity,  temperature  prohles  and  their  derivatives  have  shown 
that  they  agree  well  with  the  self-similar  solution  for  a  bonndary  layer  over  a  flat  plate 
(see  Appendix  A.l).  Only  the  second  derivatives  of  the  velocity  and  temperature  prohles 
demonstrate  some  differences  between  the  DNS  results  and  the  self-similar  solution.  Because 
of  the  viscous-inviscid  interaction,  the  edge  velocity,  Ue,  and  temperature,  T^,  are  slightly 
different  from  the  free-stream  values  Uoo  and  T^,  respectively.  Therefore,  the  local  edge 
Mach  number  also  is  different  from  the  free  stream  Mach  nnmber.  However,  the  viscous- 
inviscid  interaction  is  weak  at  the  considered  how  parameters.  For  example,  we  consider  the 
how  parameters  ontside  the  boundary  layer  at  y/L  ^  50,  where  the  length  scale  L  is  dehned 
as  {yooX / PooUacY^"^ ■  At  a;  =  0.099  m,  0.359  m,  and  0.659  m,  the  dimensionless  velocities 
and  the  local  Mach  nnmbers  are  equal  to  U/Uoo  =  0.9983,  0.9992,  0.9994,  and  =  5.845, 
5.882,  5.892,  respectively.  Therefore,  we  neglect  these  small  variations  in  the  local  Mach 
nnmber  Mg,  in  the  edge  velocity,  Ue,  and  temperature,  Tg,  and  consider  them  to  be  equal  to 
the  free-stream  parameters. 

Comparison  of  the  eigenvalues  a  obtained  using  the  self-similar  and  DNS  prohles  is  shown 
in  Appendix  A.l.  There  is  a  diherence  in  at  high  frequencies.  In  the  analysis  of  the  hat 
plate  data,  the  self-similar  prohles  have  been  nsed  in  the  stability  equations.  The  analysis 
of  the  perturbations  is  limited  to  the  DNS  data  corresponding  to  pertnrbations  of  100  kHz 
only.  The  corresponding  dimensionless  freqnency  is  55.02939  x  10“®.  In  order  to  deal  with 
the  two-dimensional  perturbations  within  the  solver  of  Refs.  [Tnm07,  Tnm08],  the  spanwise 
wave  nnmber  (3  scaled  with  the  Blasius  length  scale,  L,  was  chosen  equal  to  10“®. 

In  order  to  illnstrate  further  analysis  of  DNS  results,  features  of  the  spectrum  should 
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Figure  2.3:  Real  parts  of  the  phase  velocities  of  the  discrete  modes  F  and  S  scaled  with  the 
free-stream  velocity  Uoo- 


be  introduced.  Figure  2.2  shows  the  branches  of  the  continuous  spectrum  and  two  discrete 
modes  at  a;  =  0.08  m  (i?  =  1063.02).  One  of  the  discrete  modes  is  labeled  mode  F  (fast);  the 
other  is  labeled  mode  S  (slow).  The  mode  names  stem  from  their  phase  velocity  features  in 
the  vicinity  of  the  leading  edge.  One  can  see  in  Fig.  2.3  that  mode  S  is  synchronized  with  the 
slow  acoustic  wave  (c^  =  1  —  1/Moo),  whereas  mode  F  is  synchronized  with  the  fast  acoustic 
wave  {cr  =  1  +  1/Moo).  At  the  chosen  flow  parameters,  mode  F  is  always  stable,  and  mode 
S  is  the  unstable  mode.  One  can  see  that  mode  F  is  synchronized  with  the  vorticity/entropy 
modes  having  dimensionless  phase  velocity  =  1  at  a;  ~  0.25  m  (R  =  1879.17).  The 
signihcance  of  the  decaying  mode  F  stems  from  its  synchronization  with  mode  S,  where  the 
decaying  mode  can  give  rise  to  the  unstable  mode  (switching  of  the  modes),  which  may  lead 
to  transition  [FKOl]. 

Mode  S. 

Fignre  2.4  shows  the  pressnre  pertnrbation  on  the  wall  (scaled  with  the  free-stream 
pressnre)  obtained  in  the  DNS  and  projections  onto  the  discrete  mode  S.  Amplification  of 
the  discrete  mode  evalnated  with  and  withont  the  nonparallel  flow  effects  (MMS  and  LST, 
respectively)  is  also  presented  in  Fig.  2.4.  One  can  see  that  the  nonparallel  flow  effect  is 
signihcant  in  this  example.  The  DNS  data  for  the  wall  pressnre  pertnrbation  have  wiggles 
near  the  actnator  region  dne  to  inpnt  from  the  various  modes  presented  in  the  signal  (Fig. 
2.5).  The  hltered  out  amplitude  of  the  unstable  mode  S  is  smooth,  and  it  is  in  good  agreement 
with  the  theoretical  prediction  on  the  whole  interval. 

Mode  F 

It  is  interesting  to  look  at  the  hltered  ont  decaying  mode  F  on  Fig  2.6.  It  is  in  good 
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Figure  2.4:  Projection  of  the  DNS  results  onto  the  discrete  mode  S. 


Figure  2.5:  A  closer  view  of  the  results  in  Fig.  2.4  in  the  vicinity  of  the  actuator. 
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agreement  with  the  theoretical  prediction  up  to  a;  ~  0.25  m.  After  that,  it  experiences  a 
jump  and  the  amplitude  becomes  comparable  with  the  amplitude  of  mode  S.  The  result  can 
be  attributed  to  the  next  term  in  the  expansion  (2.6).  The  second  term,  i^iy)i  can  be 
expanded  into  the  eigenfunction  system.  It  is  the  standard  problem  of  hnding  eigenfunctions 
of  a  perturbed  operator  using  the  unperturbed  basis  [FriOO].  For  the  non- resonance  case 
when  eigenvalues  of  modes  F  and  S  are  distinct  7^  0;^),  it  is  straightforward  to  hnd  a 
projection  of  (X,  y)  onto  Aapfs  (X,  y)  (indices  S  and  F  indicate  slow  and  fast  discrete 
modes,  respectively). 

After  applying  of  the  Fourier  transform  to  the  linearized  equations  (2.1)  with  respect 
to  coordinate  2:  and  substitution  of  {x,X,y),  one  can  derive  the  following  equation  for 

ASq.Vs): 


dy 


(1)' 


dy 


(1) 


dy 


_  dDs  (X) 


dX 


H2A 


(0) 

s 


(0) 


+  A 


( Yi  TT 


+  Ds  (X)  H4A 


(0) 
s  > 


(2.19) 


where  H4  =  e“^H4.  For  the  purpose  of  brevity,  we  use  subscript  S  to  indicate  the  slow 
discrete  mode  having  wave  numbers  as-  One  can  represent  the  solution  for  A^^^  as  an 
expansion  into  the  eigenfunctions  of  the  undisturbed  operator.  In  symbolic  form,  we  write 


A?=  5:  a(X)A[°K 


(2.20) 


The  symbolic  form  of  the  expansion  (2.20)  means  that  we  include  expansion  into  the  discrete 
modes  and  continuous  spectrum  as  well.  Assuming  that  there  is  no  resonance  {a^  7^  as), 
one  can  substitute  A^^^  from  Eq.  (2.20)  into  Eq.  (2.19).  Using  the  dot  product  with  the 
adjoint  eigenvector  'Bp\  we  arrive  at  the  coefficient  Cp: 


Cf{X) 


Ds  (A) 

i  [ap  —  as) 


dX  ’ 


H4A 


(0) 
s  > 


(2.21) 


The  input  of  mode  F  into  the  second  term  of  Eq.  (2.6)  has  a  wave  number  (and  phase  speed) 
corresponding  to  mode  S.  We  refer  to  this  contribution  of  the  mode  F  as  “S2F  centaur” 
in  order  to  emphasize  the  two-fold  character  of  the  term.  The  wall  pressure  perturbation 
associated  with  S2F  centaur  is  shown  in  Fig.  2.7.  Although  the  theoretical  result  for  mode  F 
downstream  from  the  point  of  synchronism  demonstrates  qualitatively  the  same  qualitative 
behavior  as  the  amplitude  of  the  DNS  projection  onto  the  mode  F,  there  is  a  quantitative 
discrepancy  that  has  yet  to  be  understood. 

Velocity  profiles  of  modes  S  and  F 

Having  found  the  coefficients  in  the  projection  of  the  DNS  results  onto  modes  S  and  F, 
we  can  compare  the  velocity  prohles  of  the  modes  with  the  DNS  results  in  order  to  evaluate 
their  signihcance  at  different  distances  from  the  actuator. 
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Figure  2.6:  Projection  of  the  DNS  results  onto  the  discrete  mode  F. 


Figure  2.7:  Projection  of  the  DNS  results  onto  the  discrete  mode  F,  and  amplitude  of  S2F 
centaur. 
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Figure  2.8:  Streamwise  velocity  perturbation  at  a;  =  0.219  m. 


y/L 


Figure  2.9:  Streamwise  velocity  perturbation  at  a;  =  0.334  m. 
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Figure  2.10:  Streamwise  velocity  perturbation  at  a;  =  0.659  m. 


One  can  see  that  the  main  input  into  the  velocity  perturbation  stems  from  mode  S.  In 
the  vicinity  of  the  modes’  synchronism,  the  amplitude  of  mode  F  is  higher  as  expected  from 
results  presented  in  Fig.  2.7.  The  receptivity  studies  in  Ref.  [TWZ07]  showed  that  the 
amplitude  of  mode  F  in  the  vicinity  of  the  actuator  is  higher  than  the  amplitude  of  mode 
S.  In  the  present  work,  we  consider  perturbations  far  downstream  from  the  actuator,  and 
mode  S  has  a  larger  amplitude  than  mode  F. 

2.4.2  Sharp  wedge 

In  this  example,  periodic-in-time  wall  blowing-suction  was  introduced  into  the  boundary 
layer  over  a  wedge  of  half-angle  5.3  degrees.  The  free-stream  flow  conditions  were:  Mach 
number  M^o  =  8,  temperature  Too  =  54.8  K,  pressure  Poo  =  389  Pa.  The  Prandtl  number  and 
the  specihc  heats  ratio  are  0.72  and  1.4,  respectively.  The  viscosity  coefficient  is  calculated  by 
Sutherland’s  law.  The  flow  parameters  are  the  same  as  in  the  work  by  Malik  et  al.  [MLS99]. 
Results  of  Ref.  [MLS99]  were  used  for  the  validation  of  the  code  used  in  the  present  work 
(see  Ref.  [TWZ07]). 

The  periodic-in-time  blowing-suction  slot  has  coordinates  of  the  leading  and  trailing 
edges  at  x*  =  51.84  mm  and  Xg  =  63.84  mm,  respectively.  These  flow  parameters  and  the 
actuator  location  correspond  to  Case  3  considered  in  Ref.  [TWZ07].  The  dimensionless  wall 
blowing-suction  of  mass  flux  at  the  wall  is  expressed  as 

15 

(pv)'  =  qog{l)  ^  sin(a;„t),  (2.22) 

71  =  1 
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Figure  2.11:  Pressure  perturbations  (pwaii/Poo)  on  the  wedge  at  three  frequencies. 

where  go  =  0.734  x  10“^  is  a  dimensionless  amplitude  parameter,  scaled  by  the  free-stream 
streamwise  mass  flux;  function  g{l)  is  dehned  in  Eq.  (2.17),  and  ujn  is  the  circular  frequency 
of  the  multi-frequency  perturbations. 

Figure  2.11  illustrates  the  pressure  perturbations  on  the  wedge  at  three  frequencies:  44.76 
kHz,  104.44  kHz,  and  164.12  kHz  (the  dimensionless  frequencies  upe/iPeU^)  scaled  with 
the  local  boundary-layer  edge  parameters  are  approximately  21.4  x  10“®,  50.0  x  10“®,  and 
78.6  X  10“®,  respectively).  Figure  2.12  shows  the  local  Reynolds  number  R  =  {pf,x/ pJJeY^'^ 
versus  coordinate  x. 

In  the  following  examples,  analysis  of  the  flow  stability  was  based  on  the  velocity  and 
temperature  prohles  obtained  from  the  computations  without  an  assumption  about  the  self¬ 
similar  character  of  the  boundary  layer  flow.  In  order  to  compare  the  projected  amplitudes 
with  those  predicted  using  the  method  of  multiple  scales,  we  need  derivatives  d/dx  of  the 
streamwise  velocity  and  temperature  prohles  of  the  mean  how.  These  derivatives  were  derived 
using  the  computational  prohles  together  with  the  assumption  that  the  prohles  locally  are 
self-similar. 

Figure  2.13  shows  the  imaginary  part  of  the  wave  number  a  scaled  with  L  =  {pooX / PooUqoY^‘^ 
obtained  using  the  quasi-parallel  approximation  (LST)  and  using  the  method  of  multiple 
scales  (MMS)  for  perturbations  oi  F  =  44.6  kHz.  Figure  2.14  shows  wall  pressure  pertur¬ 
bations  in  DNS  results  and  in  their  projection  onto  mode  S  together  with  the  theoretical 
prediction  when  the  nonparallel  how  ehects  are  included.  Figures  2.15  and  2.16  show  similar 
results  corresponding  to  frequency  /  =  104.44  kHz.  Figs.  2.17  and  2.18  demonstrate  the 
results  at  frequency  /  =  164.12  kHz.  It  is  interesting  that  the  nonparallel  how  ehect  on 
ttjL  has  a  diherent  character  at  low  and  high  local  Reynolds  numbers  (see  Figs.  2.15  and 
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Figure  2.12:  The  local  Reynolds  number  versus  coordinate  x. 


2.17).  Using  Mack’s  terminology  [Mac69],  mode  S  is  associated  with  the  ”£rst  Mack’s  mode” 
and  the  ’’second  Mack’s  mode”  at  low  and  high  local  Reynolds  numbers,  respectively  (see 
clarification  of  the  terminology  in  Ref.  [FTIO]).  The  nonparallel  flow  effects  destabilize  the 
first  mode  and  stabilize  the  second  mode.  The  same  observation  was  reported  by  Chang  and 
Malik  [CM93].  Appearance  of  the  strong  hump  in  aiL  at  low  Reynolds  numbers  as  seen  in 
Figs.  2.15  and  2.17  is  attributed  to  the  transformation  of  the  first  mode  into  the  second  one 
[CM93]. 

2.5  Discussion  of  the  results 

The  presented  results  illustrate  how  the  multimode  decomposition  technique  may  serve  as  a 
tool  for  gaining  insight  into  the  flow  dynamics  in  the  presence  of  perturbations  belonging  to 
different  modes.  In  the  past,  one  could  compare  DNS  results  with  theoretical  prediction  for 
the  unstable  mode  only  far  downstream  from  an  actuator  where  the  unstable  mode  dominates 
the  total  signal.  Using  the  biorthogonal  eigenfunction  system,  one  can  compare  DNS  results 
with  theoretical  predictions  for  the  unstable  and  stable  modes  in  the  vicinity  of  the  actuator 
as  well. 

In  Ref.  [TWZ07]  and  in  the  present  work,  we  have  found  that  the  multimode  decompo¬ 
sition  requires  a  more  elaborate  analysis  within  the  point  of  synchronism  of  mode  F  with 
the  continuous  spectrum.  Therefore,  an  extension  of  the  theoretical  model  of  Ref.  [FKOl] 
to  the  case  of  continuous  spectrum  is  required. 

Analysis  of  the  growth  rates  using  the  method  of  multiple  scales  revealed  that  they  have 
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Figure  2.13:  Im{ai)  versus  x  at  frequency  /  =  44.76  kHz. 


f  =  44.76  kHz 


Figure  2.14:  Projection  of  the  DNS  results  onto  the  discrete  mode  S  and  comparison  with 
the  theoretical  prediction  using  the  method  of  multiple  scales  (MMS).  /  =  44.76  kHz. 
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Figure  2.15:  Im{ai)  versus  x  at  frequency  /  =  104.44  kHz. 
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a  non-monotonic  character  in  the  region  where  the  hrst  Mack’s  mode  is  transformed  into 
the  second  one  (the  terminology  regarding  Mack’s  modes  is  clarihed  in  Ref.  [FTIO]).  This 
observation  is  in  agreement  with  studies  in  Ref.  [CM93].  The  nonparallel  flow  effects  in  the 
boundary  layer  over  an  adiabatic  wall  destabilize  the  first  mode  and  stabilize  the  second  one. 
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f  =  104.44  kHz 
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Figure  2.16:  Projection  of  the  DNS  results  onto  the  discrete  mode  S  and  comparison  with 
the  theoretical  prediction  using  the  method  of  multiple  scales  (MMS).  /  =  104.44  kHz. 
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Figure  2.17:  Im{ai)  versus  x  at  frequency  /  =  164.12  kHz. 
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f  =  164.12  kHz 


Figure  2.18:  Projection  of  the  DNS  results  onto  the  discrete  mode  S  and  comparison  with 
the  theoretical  prediction  using  the  method  of  multiple  scales  (MMS).  /  =  164.12  kHz. 
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Chapter  3 


High-speed  boundary- layer  instability: 
Old  terminology  and  a  new  framework 


This  section  is  based  on  the  paper  published  in  collabaration  with  A.  Fedorov  in  Refs.  [FTll] 
and  [FTIO] 

3.1  Introduction 

In  hydrodynamic  stability  analysis,  the  normal  mode  concept  suggests  considering  sim¬ 
ple  wave-like  solutions  of  the  linearized  Navier-Stokes  equations  in  the  form  q'{x,y,t)  = 
q{y)  exp{iax  —  iut),  where  q{y)  is  a  complex  amplitude  function.  This  form  assumes  that 
the  mean  flow  is  parallel  or  quasi-parallel;  i.e.  the  mean  flow  velocity  and  temperature  pro- 
hles  are  functions  of  y  only.  The  underlying  idea  of  the  normal  mode  concept  is  that  instead 
of  solving  a  specific  physical  initial  boundary-value  problem,  we  consider  normal  modes  of 
the  wave  system,  expecting  that  a  solution  of  the  initial  boundary-value  problem  for  the 
partial  differential  equations  (PDF)  can  be  presented  as  a  sum  of  the  normal  modes.  If 
there  is  an  unstable  mode,  it  is  expected  that  this  mode  will  be  present  in  the  solution  of  a 
specihc  physical  problem  (realized  experimentally  or  solved  numerically)  and  it  can  be  dom¬ 
inant  after  sufficient  amplihcation  in  space  and/or  time.  Because  the  continuum  medium 
has  an  inhnite  number  of  degrees  of  freedom,  the  disturbance  held  has  a  numerable  discrete 
spectrum  or  a  continuous  spectrum,  or  a  combination  of  both  [CJJ67,  Mik74]. 

In  general,  we  may  consider  the  wave  number,  a  =  and  frequency,  u  =  Ur  +  icUi, 

as  complex  quantities.  Very  often  two  formulations  are  mentioned:  temporal  and  spatial 
theories.  In  the  temporal  stability  theory,  Oj  =  0  and  is  considered  as  a  parameter,  while 
u  is  the  complex  frequency  determined  from  the  dispersion  relation  u  =  u!{a,  R,  M, ...). 
According  to  the  normal  mode  treatment,  the  solution  depends  on  time  as  exp(— icut)  = 
exp{—iujrt  +  oJit).  It  is  unstable  (grows  exponentially  with  t)  if  ooi  >  0,  stable  if  uii  <  0,  and 
neutral  if  cuj  =  0.  In  the  spatial  stability  theory,  uji  =  Q  and  ojr  is  the  real  frequency  con¬ 
sidered  as  a  parameter,  while  the  complex  wavenumber  a  is  determined  from  the  dispersion 
relation  a  =  a{u,  R,  M, ...).  Now,  the  disturbance  amplitude  depends  on  the  streamwise 
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coordinate  as  exp(iaa;)  =  eyiY>{iarX  —  aix).  Propagating  downstream,  the  disturbance  grows 
exponentially  with  a;  if  o;*  <  0,  and  decays  if  a*  >  0. 

One  can  hnd  in  the  literature  an  ambiguity  with  the  choice  between  temporal  or  spatial 
formulations  of  stability  problems.  The  ambiguity  has  its  origin  at  the  step  when  partial 
solutions  of  PDE  are  suggested  in  the  form  of  normal  modes.  At  this  point,  the  link  with  the 
initial  boundary- value  problem  for  a  specihc  setup  is  lost,  and  we  have  an  artificial  dilemma 
of  spatial  versus  temporal  formulations. 

Gustavsson  [Gus79]  solved  an  initial  value  problem  for  perturbations  in  the  incompress¬ 
ible  boundary  layer.  Salwen  and  Grosch  [SG81]  proved  that  Gustavsson  solution  is  an 
expansion  into  the  modes  of  discrete  and  continuous  spectra  of  the  temporal  stability  prob¬ 
lem.  The  difference  between  the  discrete  spectrum  and  the  continuous  spectrum  is  their 
behavior  outside  the  boundary  layer  {y  cx)).  The  discrete  spectrum  is  required  to  vanish 
as  1/  — cx),  while  the  continuous  part  is  only  required  to  be  bounded.  Ashpis  and  Reshotko 
[AR90]  considered  the  signaling  problem  when  a  periodic-in-time  point  source  introduces 
perturbations  through  the  wall  starting  at  t  =  0.  They  showed  that  the  solution  of  this 
problem  can  be  presented  (after  a  sufficiently  long  time  when  the  transient  effect  is  smeared 
out)  as  an  expansion  into  modes  of  the  discrete  and  continuous  spectra  resulted  from  the 
spatial  stability  theory.  Therefore,  the  analysis  of  a  specihc  initial  boundary-value  problem 
for  PDE  leads  to  the  unique  expansion  of  a  solution  into  the  normal  modes  without  any  am¬ 
biguity.  In  compressible  boundary  layers,  there  are  seven  branches  of  continuous  spectrum 
[TF83,  BM92]:  three  branches  correspond  to  waves  propagating  upstream  with  rapid  decay, 
two  branches  correspond  to  slow  and  fast  acoustic  waves  propagating  downstream,  and  two 
branches  correspond  to  vorticity  and  entropy  waves  propagating  downstream. 

Mack  [Mac69]  discovered  new  unstable  modes  of  the  discrete  spectrum  in  high-speed 
boundary  layers.  He  hrst  considered  the  temporal  stability  problem  for  inviscid  perturba¬ 
tions  and  computed  the  eigenvalue  dependencies  uj{a)  for  the  hat-plate  boundary  layer  at 
various  free-stream  parameters.  In  discussion  of  these  dependencies.  Mack  used  the  terms 
modes  and  families.  However,  his  usage  of  the  term  modes  is  inconsistent  with  the  mathemat¬ 
ical  dehnition  of  modes  as  eigenfunctions  that  can  be  used  for  expansion  of  PDE  solutions. 
Mack  dehnition  of  families  is  what  mathematicians  would  normally  call  modes.  Moreover, 
Mack  subdivided  a  family  representing  one  dependency  of  uj{a)  into  diherent  pieces  giving 
them  diherent  labels  of  modes.  To  illustrate  Mack  terminology,  we  refer  to  Fig.  3.1  showing 
the  phase  speed  diagram  Ci{a)  =  Ui{a)/a  for  the  boundary  layer  on  a  thermally  insulated 
(adiabatic)  hat  plate  at  the  Mach  number  M  =  3.8.  One  may  compare  Fig.  3.1  with 
Fig.  11.6a  of  [Mac69].  Note  that  because  we  conduct  computations  at  constant  Prandtl 
number,  Pr  =  0.72,  and  in  some  cases  at  diherent  stagnation  temperature,  there  are  some 
diherences  between  our  results  and  the  results  of  [Mac69].  Hereafter,  dimensionless  param¬ 
eters  are  obtained  using  the  Blasius  length  scale  a/ u^x* /U*  (asterisk  denotes  a  dimensional 
quantity).  Following  Mack  notation,  the  wavenumbers  asn  in  =  1,  2, ...)  are  associated  with 
the  generalized  inhection  point  at  which  d  {{dU/dy)  /T)  /dy  =  and  the  wavenumbers  ain 
{n  =  1,  2, ...)  are  associated  with  regular  neutral  solutions. 

Mack  wrote  (page  11-21,  [Mac69]):  The  eigenvalues  for  q  7^  0  he  on  two  separate  curves 
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Figure  3.1:  Eigenvalue  diagramm.  M  =  3.8,  Tq  =  300K,  7  =  1.4,  Pr  =  0.72,  adiabatic 
wall.  Inviscid  theory.  Region  A  corresponds  to  first-mode  amplified  solutions.  Region 
B  corresponds  to  second-mode  amplified  solutions.  Region  C  corresponds  to  second-mode 
damped  solutions. 


in  Fig.  11.6  and  form  distinct  families  of  solutions.  One  family  starts  at  the  neutral  sonic 
solution;  the  other  starts  at  the  hrst-mode  regular  neutral  solutions.  It  is  easier  to  discuss  the 
inviscid  solutions  in  terms  of  these  families  rather  than  in  terms  of  the  modes.  The  mode 
identihcation  changes  along  the  eigenvalue  curves,  and  is  often  not  well  dehned.  Further 
down  the  same  page.  Mack  introduces  his  dehnitions:  In  Fig.  11.6  the  solutions  of  the  cq 
family  between  (0,  cq)  and  (a^i,  Cg)  are  amplihed,  and  there  is  almost  no  phase  change  across 
the  boundary  layer  in  the  pressure  fluctuations.  These  solutions  will  be  called  hrst-mode 
amplihed  solutions.  The  solutions  of  the  an  family  between  (an,!)  and  {as2,Cs)  are  also 
amplihed,  but  except  near  (an,  1)  there  is  a  sizable  phase  change  in  the  pressure  huctuations. 
These  solutions  will  be  called  second-mode  amplihed  solutions.  In  the  next  paragraph  on  the 
same  page.  Mack  introduces  a  dehnition  of  the  other  part  of  the  cq  family:  It  is  of  interest 
that  there  is  a  local  maximum  in  the  damping  rate  of  the  Cq  family  of  solutions  at  almost 
the  same  wave  number  for  which  the  amplihcation  rate  of  the  solution  of  the  an  family  is 
a  maximum.  Since  the  pressure-huctuation  phase  change  is  large  for  the  damped  solutions 
near  this  wave  number,  those  solutions  will  be  called  second-mode  damped  solutions.  Thus, 
Mack  uses  two  diherent  mode  labels  for  one  family  curve. 

The  spectrum  structure  becomes  more  complicated  at  higher  Mach  numbers.  For  exam¬ 
ple,  Figure  3.2  shows  Cj(a)  at  M  =  7  (it  is  similar  to  Fig.  11.8a  of  [Mac69]).  In  connection 
with  this.  Mack  wrote  on  page  11-23  of  [Mac69]:  As  shown  in  Fig.  11.8a  for  M  =  7,  the 
Co  family  of  solutions  no  longer  contains  the  hrst-mode  neutral  subsonic  solution  (asi,Cs) 
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Figure  3.2;  Eigenvalue  diagramm.  M  =  7,  Tq  =  540K,  7  =  1.4,  Pr  =  0.72,  and  Tg  =  50K. 
Region  A  corresponds  to  first-mode  amplified  solutions.  Region  B  corresponds  to  second¬ 
mode  amplified  solutions.  Region  C  corresponds  to  second-mode  damped  solutions.  Region 
D  corresponds  to  third-mode  damped  solutions. 
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and  the  second-mode  damped  solntions,  bnt  instead  contains  the  second-mode  amplihed 
solntions,  the  second-mode  neutral  subsonic  solution  {as2,Cs),  and  the  third-mode  damped 
solutions.  These  examples  illustrate  that  Mack  dehnitions  of  modes  are  inconsistent  with 
conventional  usage  of  the  term  normal  modes.  In  Figs.  3.1  and  3.2,  two  distinct  normal 
modes  correspond  to  the  cq  and  (an,  1)  families. 

Mack  then  considered  viscous  perturbations  and  compared  the  dependencies  of  c(a)  at 
finite  Reynolds  numbers  with  the  inviscid  ones.  For  one  of  Macks  choices  of  parameters, 
there  were  two  separate  inviscid  amplification  rate  curves  for  the  first  and  second  modes 
([Mac69],  page.  12-24)  (i.e.,  two  inviscid  normal  modes),  but  only  a  single  amplification  rate 
curve  at  the  hnite  Reynolds  number  shown  (i.e.,  one  viscous  normal  mode).  This  example 
illustrates  that  the  spectrum  structure  can  be  drastically  changed  by  the  viscous  effects. 

Receptivity  studies  in  1980s-2000s,  which  are  summarized  in  [Fed03,  Tum06a],  indicated 
that  excitation  of  instabilities  predominantly  occurs  upstream  from  the  unstable  region. 
This  motivated  Fedorov  and  Khokhlov  [FKOl]  to  investigate  the  disturbance  spectrum  in  the 
region  located  upstream  from  the  lower  neutral  branch.  The  spatial  stability  analysis  [FKOl] 
of  2D  unstable  modes  in  hypersonic  flows  revealed  the  following  features  (Fig.  3.3):  (1)  in  the 
leading-edge  region,  two  discrete  modes  (fast  and  slow)  are  synchronized  with  the  fast  and 
slow  acoustic  waves  of  the  continuous  spectrum,  respectively;  (2)  further  downstream,  the 
fast  discrete  mode  is  synchronized  with  the  entropy  and  vorticity  waves  of  the  continuous 
spectrum;  (3)  further  downstream,  the  fast  discrete  mode  is  synchronized  with  the  slow 
discrete  mode.  The  latter  synchronization  leads  to  a  branching  of  the  discrete  spectrum 
[GF,  GF90].  As  shown  in  [FKOl],  the  normal  mode  decomposition  is  not  valid  in  the  branch¬ 
point  vicinity  and  should  be  replaced  by  a  local  solution  accounting  for  the  coupling  of 
discrete  modes.  Gonsequences  of  these  hndings  for  modeling  of  high-speed  transition  are 
recently  discussed  by  Fedorov  [Fedll]. 

The  aforementioned  shortcomings  of  Mack  terminology  motivated  Fedorov  [Fed03]  to 
categorize  the  discrete  modes  using  their  asymptotic  behavior  near  the  leading  edge.  As 
shown  in  [FKOl]  for  2D  disturbances,  the  phase  speed  of  one  mode  tends  to  c  =  1  —  1/M 
of  slow  acoustic  wave,  whereas  the  phase  speed  of  the  other  tends  toc=l-|-l/Mof  fast 
acoustic  wave.  Fedorov  [Fed03]  called  the  former  mode  slow  mode  (mode  S)  and  the  latter 
fast  mode  (mode  F).  This  terminology  came  to  be  used  for  the  interpretation  of  DNS  results 
[EFS06,  TWZ07,  WZ23,  EFS08]  on  the  receptivity  and  stability  of  high-speed  boundary- 
layer  flows. 

The  temporal  stability  analyses  [FA03,  FT05]  revealed  similar  features  of  the  spectral 
curves  u{a)  at  hxed  Reynolds  number  R.  Namely,  the  fast  (or  slow)  mode  is  synchronized 
with  the  fast  (or  slow)  acoustic  wave  in  the  long-wave  limit  a  — 0.  As  a  increases,  the  fast 
mode  is  synchronized  with  the  entropy /vorticity  waves,  and  then  it  is  synchronized  with  the 
slow  mode  that  leads  to  branching  of  the  discrete  spectrum.  The  branching  causes  instability 
of  what  is  called  the  second  mode  (according  to  Mack  terminology).  This  similarity  is  not 
surprising  because  the  spectrum  behavior  depends  on  the  product  a*6.  The  case  of  a*  =£xed 
and  (5  — 0  is  treated  as  a  leading-edge  limit,  while  the  case  of  6  =fixed  and  a*  — Ois  treated 
as  a  long- wave  limit. 


Mack  [Mac69]  introduced  the  viscous  and  inviscid  families  of  solutions  when  the  recep¬ 
tivity  problem  was  not  understood,  and  the  decomposition  of  the  solutions  of  the  linearized 
Navier-Stokes  equations  had  not  been  developed.  A  loose  usage  of  the  hst  mode  and  the 
second  mode  terminology  has  not  been  an  obstacle  in  discussion  of  LST  applications  to  tran¬ 
sition  prediction  (e.g.  using  the  e^-method)  even  in  the  cases  when  only  one  unstable  normal 
mode  exists  at  hnite  Reynolds  numbers.  However,  one  has  to  be  aware  of  the  meaning  of 
this  imprecise  terminology  when  transition  prediction  is  based  on  consideration  of  the  hrst 
and/or  second  mode.  On  the  other  hand,  the  usage  of  the  mode  F  and  mode  S  terminol¬ 
ogy  has  not  yet  been  discussed  properly  in  connection  with  the  inviscid  stability  theory.  In 
this  paper,  we  consider  the  discrete  spectrum  of  disturbances  in  high-speed  boundary  layers 
in  order  to  clarify  the  contemporary  terminologies  and  their  interrelations.  In  addition,  a 
model  of  discrete  spectrum  branching  is  outlined  in  order  to  explain  different  behaviors  of 
the  aforementioned  discrete  modes. 


3.2  Discussion  of  discrete  spectrum  in  high-speed 
boundary  layers 

This  section  illustrates  examples  of  discrete  normal  modes  resulted  from  numerical  solutions 
of  the  temporal  and  spatial  stability  problems  for  the  locally  parallel  boundary- layer  flow  on 
a  flat  plate.  All  computations  were  conducted  for  perfect  gas  of  constant  Prandtl  number 
Pr  =  0.72  and  specihc  heat  ratio  7  =  1.4.  Viscosity  is  calculated  using  the  Sutherland  law 
with  a  constant  of  110. 4K  and  the  bulk  viscosity  is  zero.  The  mean  flow  prohles  correspond 
to  the  compressible  Blasius  solution;  i.e.,  the  viscos-inviscid  interaction  is  neglected.  The 
discrete  spectrum  structures  for  viscous  (at  hnite  Reynolds  number)  and  inviscid  {R  00) 
disturbances  are  compared. 

3.2.1  Temporal  stability  theory 

In  order  to  show  the  link  between  discrete  modes  at  hnite  and  inhnite  Reynolds  numbers, 
we  have  to  consider  the  regular  neutral  modes.  The  inviscid  stability  equations  have  neutral 
(cj  =  0)  discrete  solutions  with  the  phase  speed  1  <  <  1  -|-  1/M.  Figures  3.4a  and  3.4b 

illustrate  Ci{a)  and  Cr{a)  for  the  boundary  layer  on  an  adiabatic  hat  plate  at  M  =  3.8  and 
R  ^  00,  while  hgures  3.5a  and  3.5b  show  Ci{a)  and  Cr(a)  at  hnite  R  =  2000.  In  the  viscous 
case,  the  fast  discrete  mode,  which  is  synchronized  with  fast  acoustic  wave  of  c  =  1  -|-  1/M 
as  a  — 0,  has  a  small  jump  of  Ci{a)  at  the  point  where  Cr  =  1.  Actually,  we  should  treat 
the  fast  mode  as  two  discrete  modes.  One  of  them  has  >  1  and  coalescences  with  the 
continuous  spectrum  branch  corresponding  to  vorticity  and  entropy  waves  (c^  — 1  -|-  0). 

The  other  mode  departs  from  this  branch  with  Cr  1  —  0  and  has  Cr  <  1.  The  coalescence 
and  departure  occur  at  slightly  diherent  values  of  q.  Because  the  jump  of  q  is  small,  it  is 
usually  ignored  in  stability  analyses.  However,  the  eigenfunctions  are  substantially  different 
on  different  sides  of  =  1.  Their  outer  asymptotic  forms  contain  different  eigenvectors, 
which  are  proportional  to  exp(iax  ±  i\y  —  iuit)  as  1  =F  0  and  y  ^  00.  This  difference 
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b) 

Figure  3.4:  c*  (a)  and  (b)  versus  a.  M  =  3.8,  Tq  =  300K.  Inviscid  theory. 


b) 

Figure  3.6:  c*  (a)  and  (b)  versus  a.  M  =  4.2,  Tq  =  300K.  Inviscid  theory- 
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is  important  for  analysis  of  the  initial  value  problem  [FA03]  describing  receptivity  to  vor- 
ticity/entropy  disturbances.  Hereafter  we  call  these  modes  as  mode  F-  and  mode  F+.  As 
shown  in  hgure  3.5,  mode  F-  is  associated  with  the  regular  neutral  solutions  in  the  inviscid 
limit.  Mode  S,  which  is  synchronized  with  the  slow  acoustic  wave  of  phase  speed  c  =  1  —  1/M 
as  a  — 0,  corresponds  to  the  inviscid  slow  mode  shown  in  Figs.  3.4a  and  3.4b.  Mode  F+  is 
unstable  at  the  range  of  relatively  large  relevant  to  the  second-mode  instability.  Mode  S  is 
unstable  in  the  range  of  small  relevant  to  the  first-mode  instability. 

Figures  3.6a  and  3.6b  illustrate  Ci{a)  and  Cr{a)  for  the  boundary  layer  on  an  adiabatic  flat 
plate  at  M  =  4.2  in  the  inviscid  limit  R  oo.  The  spectrum  behavior  is  qualitatively  the 
same  as  in  the  case  of  M  =  3.8  (Figs.  3.4a  and  3.4b).  However,  it  is  drastically  different  at 
the  hnite  Reynolds  number  R  =  2000  (Figs.  3.7a  and  3.7b).  Now,  there  is  only  one  unstable 
mode  (mode  S)  that  has  two  maxima  of  c*  associated  with  the  two  unstable  regions.  This 
is  similar  to  Mack  observation  in  the  cases  of  M  =  5.8  and  M  =  7  ([Mac69],  page  12-24). 
Comparing  the  cases  shown  in  Figs  3.5a,b  and  Figs.  3.5a,b  we  conclude  that  increasing  of 
Mach  number  from  3.8  to  4.2  (at  hxed  R  =  2000)  leads  to  a  qualitative  change  in  the  discrete 
spectrum  pattern. 

Figures  3.8a  and  3.8b  shows  q  and  Cr  versus  a  for  the  inviscid  and  viscous  cases  at  M  =  7. 
Comparison  of  these  data  with  that  shown  in  Figs.  3.6a  and  3.6b  indicates  that  in  the  limit 
of  i?  — cx)  the  discrete  spectrum  pattern  switches  over  to  a  different  topology  somewhere  in 
the  range  4.2  <  M  <  7.  This  change  is  similar  to  that  observed  for  viscous  disturbances  in 
the  Mach  number  range  from  3.8  to  4.2. 

3.2.2  Spatial  stability  theory 

Adiabatic  wall 


In  the  case  of  an  adiabatic  wall  in  the  spatial  theory,  available  numerical  results  indicate 
that  there  is  only  one  unstable  mode  (mode  S).  Figures  3.9a  and  3.9b  show  ai  versus  fre¬ 
quency  at  M  =  3.8  and  hxed  Reynolds  numbers.  One  can  see  that  there  are  two  minima 
in  ai{u).  At  R  =  600,  only  the  hrst  minimum  represents  spatially  unstable  solutions.  At 
R  =  2000,  unstable  perturbations  can  be  found  in  two  bands  of  frequencies.  Plotting  the 
neutral  curves  in  the  uj  —  R  plane  shows  two  unstable  regions,  but  both  of  them  are  associated 
with  mode  S. 

Another  example  is  shown  in  Fig.  3.10  for  the  case  of  M  =  4.2  and  R  =  2000.  Again, 
mode  S  is  unstable  in  the  low-frequency  band  relevant  to  the  hrst-mode  instability,  and  in 
the  high-frequency  band  relevant  to  the  second-mode  instability. 

Cold  wall 


It  was  shown  in  [FKOl]  that  mode  F-|-  can  be  unstable  as  well.  In  order  to  illustrate  the 
spectral  pattern,  we  consider  perturbations  at  a  hxed  frequency  parameter  F  = 

The  dimensionless  angular  frequency  is  cu  =  F  x  R,  the  eigenvalues  and  phase  speeds  are 
expressed  as  a  =  q;(R),  c(R)  =  u}/a{R). 


36 


b) 

Figure  3.10:  a*  (a)  and  Cr  (b)  versus  to.  M  =  4.2,  Tq  =  300K,  R  =  2000.  Adiabatic  wall. 


Figure  3.11:  ai  (a)  and  Cr  (b)  versus  R.  M  =  5.5,  Tg  =  70K,  T^/Tc 


I 


I 


Figures  3.11  and  3.12  show  ai{R)  and  Cr{R)  of  modes  S,  F-  and  F+  for  the  boundary-layer 
flow  at  M  =  5.5  and  M  =  6.5,  respectively.  In  the  both  cases,  the  mean-flow  temperature 
at  the  upper  boundary-layer  edge  is  T*  =  70K,  the  wall  temperature  is  =  Q.fT^^,  and 
the  frequency  parameter  is  F  =  10“^.  Similar  to  the  temporal  stability  cases  discussed  in 
Section  A,  mode  F-  (or  mode  S)  is  synchronized  with  the  fast  (or  slow)  acoustic  wave  in 
the  leading-edge  region  corresponding  to  the  long-wave  limit.  Then  mode  F-  is  synchronized 
with  vorticity/entropy  waves  of  the  continuous  spectrum  at  the  station  where  =  1.  As 
the  fast-mode  eigenvalue  crosses  the  vorticity/entropy  branch,  its  imaginary  part  jumps  by 
a  small  value  Aoj.  Further  downstream  mode  F_|_  is  synchronized  with  mode  S  and  the 
dispersion  curves  branch  out. 

In  the  case  of  M  =  5.5  (Figs.  3.11a  and  3.11b),  the  branching  leads  to  destabilization  of 
mode  F+  and  stabilization  of  mode  S  in  the  range  of  relevant  to  instability  of  Mack  second 
mode.  As  contrasted  to  the  adiabatic  wall  cases,  there  is  no  instability  of  the  slow  mode  at 
relatively  small  because  this  instability  (associated  with  Mack  hrst  mode)  is  suppressed  by  the 
wall  cooling.  In  the  case  of  M  =  6.5  (Figs.  3.12a  and  3.12b),  the  branching  is  qualitatively 
different.  Now  mode  S  is  unstable  and  mode  F+  is  stable.  Detailed  computations  showed 
that  the  switching  from  one  topology  to  the  other  occurs  at  M  6. 

3.3  A  model  of  the  discrete  spectrum  branching 

Guschin  and  Fedorov  [GF,  GF90]  have  realized  that  the  aforesaid  branching  of  discrete 
spectrum  is  similar  to  the  branching  of  weakly  coupled  waves  that  is  typical  for  wave  systems 
in  the  physics  of  plasmas  Fedorchenko  (see  also  64  of  [LP81]).  Herein  we  outline  an  analytical 
model  describing  this  branching. 

Gonsider  the  temporal  spectrum  uj{a)  and  assume  that  the  eigenvalues  0Ji{a)  and  U2(a) 
(in  our  case,  they  are  relevant  to  modes  F+  and  S)  are  synchronized  in  the  vicinity  of  a  certain 
point  ao'.  oji^ao)  ~  cc;2(ao)  ~  ^o-  The  other  modes  of  the  discrete  and  continuous  spectra  are 
assumed  to  be  sufficiently  far  from  the  synchronization  region.  Then,  the  dispersion  relation 
in  the  vicinity  of  (q;o,cc;o)  can  be  approximated  as 

(a)  —  aid)  {u  —  a2d)  =  St,  (3.1) 

where  u  =  u  —  uq  and  d  =  a  —  ao  are  local  variables,  et  is  a  small  parameter  characterizing 
weak  coupling  between  the  modes,  and  ai  >  02  >  0  are  the  group  velocities  of  the  decoupled 
modes  (at  et  =  0)  schematically  shown  in  Fig.  3.13.  The  dispersion  curves  branch  out  in 
the  vicinity  of  the  synchronization  point  as 


1^1,2 


_  , 

a± 


Q;2(ai  -  02)^  +  Et- 


(3.2) 


The  branching  pattern  depends  on  the  locations  of  the  branch  points  df,^i^2  =  (cq  — 

02)  in  the  complex  a-plane.  If  et  >  0,  then  d^^i  2  are  pure  imaginary  and  the  dispersion  curves 
branch  out  as  shown  schematically  in  Fig.  3.14.  In  this  case,  the  coupling  does  not  cause 
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co  =  a^a 


(o  =  co-coQ^a  =  a-a^ 


Figure  3.13:  Synchronization  of  two  uncoupled  modes  in  the  vicinity  of  (q;o,i^o)- 
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Figure  3.14:  Neutral  branching  at  St  >  0. 
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instability.  If  e*  <  0,  the  branch  points  are  real  while  chi, 2  are  complex  values  in  the  range 
abri  <  a  <  ahr2-  In  this  case,  one  mode  becomes  stable  and  the  other  unstable.  Depend¬ 
ing  on  bypasses  of  the  branch  points,  there  are  four  patterns  of  dispersion  curves  as  shown 
schematically  in  Fig.  3.15. 

For  the  temporal  spectra  of  the  boundary-layer  flow  discussed  in  Section  IIA,  the  branch¬ 
ing  of  modes  F_|_  and  S  is  topologically  identical  with  one  of  the  patterns  shown  in  Fig.  3.15. 
The  branch  points  were  determined  numerically  from  the  condition  \duj/da\~^  =  0.  Com¬ 
putations  showed  that  aferi,2  are  slightly  shifted  from  the  real  axis  to  the  lower  or  upper 
half-plane  of  complex  a  depending  on  the  basic  flow  parameters.  As  one  marches  along  real 
a,  the  branch  points  are  bypassed  from  above  or  below,  leading  to  a  particular  branching 
pattern. 

The  branch  points  move  in  the  a-plane  with  variations  of  the  flow  parameters,  say  the 
Mach  number.  As  soon  as  abr  crosses  the  real  axis,  the  branching  pattern  switches  over  from 
one  to  another.  For  example.  Figure  3.16  shows  such  switching  from  the  branching  A  to  the 
branching  B  due  to  a  slight  increase  of  Mach  number  from  4  to  4.2. 

By  computing  the  trajectories  of  the  branch  points  and  determining  the  parameters  at 
which  abri,2  crosses  the  real  axis  one  can  identify  domains  relevant  to  different  branching 
patterns.  As  an  example.  Figure  3.17  shows  the  boundary  between  the  branching  A  and 
the  branching  B  in  the  Mach  number-wall  temperature  plane.  Below  this  boundary  (at  the 
relatively  cold  wall  condition)  there  are  two  unstable  modes:  mode  S  is  unstable  in  the  range 
of  small  a  and  mode  F+  is  unstable  in  the  range  of  large  a.  Above  this  boundary,  there  is 
only  one  unstable  mode  the  slow  mode. 

Similar  analysis  can  be  conducted  for  the  spatial  stability  problem  [FKOl].  The  dispersion 
relation  is  analytically  continued  to  the  complex  plane  of  .  In  the  vicinity  of  synchronization 
point  (i?o,tto);  Equation  (1)  is  written  in  the  form 

(d  —  biu)  (d  —  b2Uj)  =  Eg,  (3.3) 

where  a  =  a  —  ao,  u  =  F{R  —  Rq)  =  FR,  bi  =  and  62  =  are  constants,  61  >  62  >  0, 
a  small  parameter  Eg  =  £4/(0102)  characterizes  the  coupling  between  synchronized  modes. 
The  solution  has  two  branches 


di,2  =  ^-^^FR  T  ^J^-{b,-b2rFm^+Eg.  (3.4) 

The  branch  points  Rbri,2  =  [F{bi  —  62)]  are  located  in  different  half-planes  of 

complex  R  depending  on  the  basic  flow  parameters.  With  marching  downstream  along  real 
R,  the  branch  points  are  bypassed  from  above  or  below  that  leads  to  different  branching 
patterns  of  a{R).  For  instance,  in  the  case  of  M  =  5.5  (Fig.  3.11),  the  both  branch  points 
he  below  the  real  axis  of  complex  i?-plane  that  gives  the  branching  B  shown  in  Fig.  3.15. 
In  the  case  of  M  =  6.5  (Fig.  3.12),  the  hrst  branch  point  is  above  the  real  axis  while  the 
second  is  below.  This  gives  the  branching  D. 

Thus,  the  analytical  model  of  weakly  coupled  synchronized  modes  helps  to  interpret 
intricate  behaviors  of  disturbance  spectrums  that  are  observed  in  high-speed  boundary  layer 
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Figure  3.15:  Branching  patterns  for  two  weakly  coupled  modes. 


a) 


b) 

Figure  3.16:  Different  branching  patterns  at  small  variation  of  Mach  number  at  Tq  =  300K, 
R  =  2000.  Branch  points  are  shown  by  black  circles  in  the  complex  a-plane.  (a)  M  =  4, 
branch  points  are  bypassed  from  below  -  branching  A,  (b)M  =  4.2,  branch  points  are 
bypassed  above  -  branching  B. 
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Figure  3.17:  ai  (a)  and  Cr  (b)  versus  R.  M  =  6.5,  Tg  =  70K,  T^/Tad  =  0.1,  F  =  10 

flows.  This  model  clearly  shows  that  the  second-mode  instability  is  due  to  synchronization 
of  modes  S  and  F+  .  Similar  mechanism  is  responsible  for  the  third-  and  higher-mode 
instabilities  in  hypersonic  boundary  layers  [GF90],  as  well  as  short-wave  instabilities  in  a 
thin  shock  layer  [GF89].  As  an  example,  figures  3.18  and  3.19  show  the  discrete  spectrum 
at  M  =  7,  Tg  =  50  K  and  R  =  20000  in  the  case  of  adiabatic  wall  (Tw  =  Tad)  and  cold 
wall  (Tu,  =  1.2Tg),  respectively.  In  both  cases,  there  are  fast  modes  f2\  F±  \  ■  ■  ■,  and  the 
slow  mode  S.  On  an  adiabatic  wall,  the  Macks  second,  third  etc.  instabilities  correspond 
to  mode  S.  These  instabilities  are  associated  with  the  branching  B  between  modes  F^'  and 
mode  S.  In  the  cold-wall  case,  the  Mack  second,  third  etc.  instabilities  correspond  to  modes 
F+\  F^^  etc.  These  instabilities  are  associated  with  the  branching  A.  As  the  number  k  of 
mode  F^^  increases,  the  growth  rates  quickly  decreases.  E.g.  the  mode  ^|3) 

in  figure  3.19a 

is  almost  neutral.  Ma  and  Zhong  [MZ03a,  MZ03b]  showed  that  the  fast  modes  f)}\  f)^^ 
etc.  (named  mode  I,  II,  etc  in  [MZ03b])  are  effectively  excited  by  fast  acoustic  waves. 

Fedorov  and  Khokhlov  [FKOl]  showed  that  the  normal  mode  decomposition  is  not  valid 
near  the  branch  point  of  discrete  spectrum.  Using  asymptotic  methods  they  obtained  a  local 
solution  providing  coupling  between  input  and  output  amplitudes  of  the  branching  modes. 
It  turned  out  that  this  coupling  is  governed  by  the  intermodal  exchange  rule,  which,  in 
the  leading-order  approximation,  does  not  depend  on  the  mean-flow  profiles.  However,  this 
result  may  be  not  valid  when  the  branch  point  is  too  close  to  the  branch  cut  of  continuous 
spectrum.  E.g.  figure  3.16  shows  that  the  relative  distance  between  the  first  branch  point 
ttbi  =  0.2199  —  iO. 00075  and  the  vorticity/entropy  branch  cut  is  Aa/a  ~  5%.  In  this  case 
the  vorticity/entropy  waves  may  significantly  affect  the  intermodal  exchange.  This  issue 
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a) 


b) 

Figure  3.18:  q  (a)  and  Cr  (b)  versus  a.  M  =  7,  =  50K,  R  =  2000.  Adiabatic  wall. 
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Figure  3.20:  Mean-flow  profiles  at  M  =  4.2,  Tq  =  300K  and  T^,  =  Tad- 


should  be  clarified  by  further  theoretical  studies. 

3.4  Eigenfunctions  of  slow  and  fast  modes 

Now  we  consider  the  eigenfunctions  of  the  slow  and  fast  modes  for  the  boundary-layer  flow 
at  M  =  4.2,  R  =  2000  and  =  Tad  (see  Figs.  3.7a,b).  In  all  cases  discussed  herein,  the 
eigenfunctions  are  normalized  to  the  absolute  value  of  the  pressure  amplitude  on  the  wall. 
Perturbations  of  pressure  p,  streamwise  velocity  u  and  temperature  6  are  made  nondimen- 
sional  using  U*  and  T*,  respectively.  Figure  3.20  shows  the  mean-flow  prohles  U{y) 

and  T{y),  where  the  boundary-layer  thickness  is  <5  =  12.5  {U (5)  =  0.99)  and  the  displacement 
thickness  is  (5i  =  9.34. 

At  q;  =  0.05  relevant  to  the  hrst  instability  maximum,  the  eigenfunction  of  mode  F_ 
typihes  regular  disturbances  of  the  phase  speed  Cr  >  1  (thin  lines  in  Fig.  3.21).  The  slow 
mode  (thick  lines)  of  Cr  =  0.81  has  the  critical  layer  in  the  vicinity  of  yc  ~  10.  Here  the 
-velocity  distribution  has  a  phase  jump  and  the  temperature  amplitude  \9\  has  a  strong 
maximum.  The  pressure  amplitude  \p\  is  almost  constant  across  the  boundary  layer  as 
typical  for  long- wave  disturbances.  Near  the  wall  (0  <  y  <  1.2)  there  is  a  well  distinguished 
viscous  sublayer  (Stokes  layer)  for  both  fast  and  slow  modes. 

Figure  3.22  shows  eigenfunctions  of  modes  F_  and  F+,  which  are  synchronized  with  the 
entropy  and  vorticity  waves  at  c,.  — >  1  ±  0.  These  functions  are  almost  identical  inside  the 
boundary  layer,  while  the  tails  for  the  u-velocity  and  temperature  are  essentially  different 
outside  the  boundary  layer.  This  clearly  shows  that  mode  F_  with  >  1  and  mode  F+  with 


50 


25 

20 

15 

10 

5 

0 

0.0  0.2  0.4  0.6  0.8  1.0  0  5  10  15  0  50  100  150  200 

iPl  |u|  lei 


Figure  3.21;  Eigenfunctions  of  the  fast  and  slow  modes  at  a  =  0.05,  M  =  4.2,  Tq  =  300K, 

Tw  =  Tad,  R  =  2000. 

Cr  <  1  are  different  normal  modes  of  the  discrete  spectrnm.  The  oscillatory  behavior  of  the 
onter  tails  is  due  to  synchronization  of  the  discrete  modes  with  the  entropy  and  vorticity 
waves,  which  oscillate  in  the  onter  flow.  A  detailed  discussion  of  these  spectrnm  featnres  is 
given  in  [FA03]. 

Figure  3.23  shows  eigenfunctions  of  modes  S  and  F+  in  the  vicinity  of  the  hrst  branch 
point  ai,ri  =  0.2199  —  f  0.00075  (see  Branching  B  in  Fig.  3.16).  As  expected,  in  this  region 
both  modes  have  similar  distribntions  across  the  boundary  layer.  However,  there  are  qnan- 
titative  differences  because  the  eigenfunctions  were  compnted  at  a  =  0.219  which  is  slightly 
shifted  from  the  branch  point. 

Fignre  3.24  illustrates  eigenfunctions  at  a  =  0.25  relevant  to  the  second  maximum  of 
instability  (see  Fig.  3.7a).  This  point  lies  approximately  in  the  middle  of  the  branching 
region  atri  <  a  <  a}jr2,  where  the  conpling  between  modes  F+  and  S  prodnces  the  most 
profonnd  effect  on  the  growth  rates.  Since  the  real  freqnencies  of  these  modes  are  practically 
identical,  their  eigenfnnctions  are  very  close  everywhere  besides  the  critical  layer,  which  is 
formed  at  the  distance  Hc  ~  10. 

The  foregoing  examples  show  that  the  eigenfnnctions  of  fast  and  slow  modes  are  qnite 
sensitive  to  singularities  of  the  dispersion  relation  caused  by  spectrum  branching  and  by 
coalescence  with  the  continuous  spectrum.  This  may  signihcantly  affect  receptivity  as  well 
as  the  nonlinear  interaction  between  the  modes. 
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Figure  3.22:  Eigenfunctions  of  fast  modes  in  the  vicinity  of  vorticity  and  entropy  waves, 
M  =  4.2,  To  =  300K,  =  Tad-,  R  =  2000.  Symbols  -  mode  T+  at  a  =  0.1899,  Cr  ^  1  —  0] 

black  lines  -  mode  T_  at  a  =  0.1886,  c,.  — 1  +  0. 
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Figure  3.23:  Eigenfunctions  of  modes  T+  and  S  in  the  vicinity  of  the  first  branch  point, 
M  =  4.2,  To  =  300K,  T^  =  Tad,  R  =  2000,  a  =  0.219. 
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Figure  3.24:  Eigenfunctions  of  fast  and  slow  modes  at  a  =  0.25,  M  =  4.2,  Tq  =  300K, 
=  Tad,  R  =  2000. 

3.5  Conclusions 

The  structure  of  the  discrete  spectrum  in  stability  analysis  of  high-speed  boundary  layers 
depends  on  basic  flow  parameters  such  as  the  Mach  number,  Prandtl  number,  Reynolds 
number  etc.  The  spectrum  may  have  two  unstable  modes  that  are  easily  associated  with 
inviscid  instabilities  of  Macks  first  and  second  modes.  With  another  choice  of  the  basic 
parameters,  the  spectrum  may  have  only  one  unstable  mode  having  two  maxima  of  the 
growth  rate. 

Terminology  introduced  by  Mack  [Mac69]  using  modes  1  and  2  (the  well-known  Erst 
and  second  modes)  is  inconsistent  with  a  proper  mathematical  treatment  of  normal  modes. 
Macks  family  of  solutions  is  actually  what  should  be  called  a  normal  mode,  whereas  his 
definitions  of  first-mode  amplihed  solutions  and  second-mode  amplihed/damped  solutions 
lead  to  confusion  when  normal  modes  are  used  to  represent  a  solution  of  the  linearized 
Navier-Stokes  equations.  The  advances  of  the  last  decades  have  made  necessary  to  develop 
new  terminology  to  replace  the  old.  In  mathematical  analysis  of  the  initial  value  problems, 
it  is  more  convenient  to  follow  the  terminology:  fast  (F)  and  slow  (S)  modes.  According 
to  this  terminology  discrete  modes  F^^  {k  =  1,2,...)  are  synchronized  with  fast  acoustic 
waves  of  the  continuous  spectrum  at  sufficiently  small  R  in  the  spatial  framework  or  at 
sufficiently  small  a  in  the  temporal  framework.  These  modes  have  >  1  and  correspond 
to  the  regular  neutral  solutions  in  the  inviscid  limit.  The  discrete  mode  S  is  synchronized 
with  slow  acoustic  wave  of  the  continuous  spectrum  as  i?  — 0  in  the  spatial  framework  or  as 
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q;  — 0  in  the  temporal  framework.  As  R  (or  a)  increases,  the  phase  speed  of  decreases, 
and  ultimately  these  modes  are  absorbed  by  the  continuous  spectrum  of  vorticity/entropy 
waves  at  Cr  =  1  +  0.  Simultaneously  modes  spring  from  the  continuous  spectrum  at 
Cr  =  1  —  0,  with  Cj  being  slightly  different  from  q  of  corresponding  modes  F^^.  Although  the 
jump  of  Cj  is  very  small  and  is  usually  ignored  in  stability  computations,  the  eigenfucntion 

(k)  ik) 

of  mode  F_^  at  =  1  +  0  is  essentially  different  from  the  eigenfucntion  of  Fy^  mode  at 
Cr  =  1  —  0.  This  difference  should  be  taken  into  account  in  solving  the  initial  value  problem 
associated  with  receptivity  to  vorticity/entropy  waves. 

Depending  on  the  flow  parameters,  both  mode  S  and  modes  F^^  can  be  unstable  (e.g.  Fig. 
19a)  or  only  mode  S  is  unstable  having  several  maxima  of  the  growth  rate  (e.g.  Fig.  18a). 
The  low-frequency  (long-wavelength)  maximum  is  associated  with  Macks  hrst  mode,  whereas 
the  higher-frequency  (shorter-wavelength)  maxima  are  associated  with  Macks  second  mode, 
third  mode  etc.  The  second- mode  instability  is  due  to  the  branching  of  modes  S  and  f!^^ 
in  the  vicinity  of  their  synchronization,  the  third-mode  instability  is  due  to  the  branching 
of  modes  S  and  Fy'^  etc.  Because  this  mechanism  quickly  weakens  as  the  mode  number 
increases,  in  many  cases  it  is  sufficient  to  consider  the  second-mode  instability  only.  For 
weakly  coupled  modes,  the  local  dispersion  relation  is  expressed  in  a  simple  analytical  form 
that  allows  us  to  explain  the  intricate  branching  patterns  observed  in  numerical  solutions. 
These  patterns  depend  on  the  locations  of  branch  points  in  the  complex  plane  (a-plane  for 
the  temporal  problem  or  i?-plane  for  the  spatial  stability  problem  describing  propagation  of 
disturbances  in  a  weakly  nonparallel  boundary  layer).  Inhnitesimal  variations  of  a  basic  flow 
parameter  (say,  Mach  number)  can  cause  a  jump  from  one  branching  pattern  to  another. 
Furthermore,  in  the  branch-point  vicinity,  the  modes  are  singular,  e.g.  their  group  velocities 
tend  to  inhnity.  Such  unphysical  behaviors  indicate  that  instead  of  isolated  normal  modes 
one  should  consider  an  initial  boundary-value  problem  formulated  for  a  certain  physical 
setup. 

To  avoid  confusion  in  terminology  we  recommend  that  authors  use  the  terms  fast  modes 
(modes  F_  and  F+)  and  slow  mode  (mode  S)  for  both  mathematical  analyses  and  the  in¬ 
terpretation  of  direct  numerical  simulations  dealing  with  initial-value  problems.  For  inter¬ 
pretation  of  experimental  data  and  LST  results  focused  on  the  characteristics  of  unstable 
disturbances,  it  is  reasonable  to  continue  to  use  Macks  terminology  with  the  terms  Macks 
first  mode  and  Macks  second  mode.  If  both  terminologies  are  involved,  the  correspondence 
between  fast  (or  slow)  modes  and  Mack  first  (or  second)  modes  should  be  established. 

The  aforementioned  singularities  of  the  discrete  spectrum  may  cause  difficulties  in  stabil¬ 
ity  analyses  based  on  traditional  LST  and  PSE  methods.  Multiple-mode  considerations  are 
needed  to  clarify  this  issue.  This  could  be  done  using  approximate  theoretical  models  such 
as  those  developed  in  [FKOl],  or  by  solving  the  linearized  Navier-Stokes  equations.  Direct 
numerical  simulations  are  also  needed  to  validate  predictions  of  the  reduced-order  models. 
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Chapter  4 

Stability  of  boundary  layers  in  binary 
mixtnres  of  oxygen  and  nitrogen 


This  section  is  based  on  the  paper  published  in  collabaration  with  J.  Klentzman  and  E.  Ulker 
[UKTll] 

NOMENCLATURE  FOR  CHAPTER  4 

c  =  complex  phase  velocity,  uj  j  a  =  Cy  ici 

Cs  =  mass  fraction  of  species  s 

Cpf  =  frozen  specihc  heat,  cal / g-mole-K 

Cp^s  =  specihc  heat  at  constant  pressure  of  species  s,  cal/g-mole 

Dij  =  binary  diffusion  coefficient,  cm? / sec 

Dij  =  pDij,  cm?  atm/ sec 

h  =  enthalpy  of  mixture,  cal /g-mole-K 

hg  =  enthalpy  of  species  s,  cal / g-mole-K 

{AIiDt^^p  =  standard  heat  of  formation  of  species  s  at  temperature  Ty^f  ,eal / g-mole 
I  =  total  enthalpy,  h  +  u^/2 

Kf  =  frozen  thermal  conductivity  of  mixture  in  thermodynamic  equilibrium,  Kf  =  Kty-\- 
Kint,  cal /cm- sec- K 

Kint  =  internal  component  of  frozen  thermal  conductivity  of  a  mixture  in  thermodynamic 
equilibrium, 

Kty  =  translational  thermal  conductivity  of  mixture  from  hrst  Chapman- Enskog  approx¬ 
imation,  cal  /  cm-sec-K 

k  =  Boltzmann’s  constant 

=  forward  reaction-rate  coefficient  for  reaction  r,  cm? / mole^-sec 
k^p  =  backward  reaction-rate  coefficient  for  reaction  r,  cm^  /mole-sec 
Le  =  frozen  binary  Lewis  number,  pCpfDij/Kf 
M  =  Mach  number 

Mg  =  molecular  weight  of  species  s,  g/g  —  mole 
Na  =  Avogadro’s  number,  molecules / g  —  mole 
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p  =  pressure,  atm 

Pr  =  frozen  Prandtl  number,  pCpf/Kf 

R  =  the  mixture  gas  constant 

Rs  =  Runiv/Mg,  gas  constant  of  species  s 

Runiv  =  universal  gas  constant,  cal/g  —  mole-K 

Sm  =  Schmidt  number,  p/pDu 

T  =  temperature,  K 

t  =  time,  sec 

u  =  x-component  of  the  velocity 
V  =  ^/-component  of  the  velocity 
w  =  ;2-component  of  the  velocity 
Wg  =  net  mass  rate  of  production  of  species  s 
X  =  coordinate  along  the  surface 
Xg  =  molar  fraction  of  species  s 
y  =  distance  from  the  wall 
=  spanwise  coordinate 
a  =  x-component  of  the  wave  number 
(3  =  ;2-component  of  the  wave  number 
5  =  boundary  layer  length  scale 

=  dehned  by  Eq.  (4.4),  cm-sec 

=  dehned  by  Eq.  (4.2),  cm-sec 

fi  =  viscosity  of  mixture  from  hrst  Chapman-Enskog  approximation,  g/ cm-sec 
p  =  density  of  mixture,  g/cm^ 

Pg  =  density  of  species  s,  g/cm^ 
u  =  complex  frequency  of  a  perturbation 

=  average  collision  cross  section  for  collisions  between  species  i  and  j,  angstrom‘s 

=  average  collision  cross  section  for  collisions  between  species  i  and  j,  angstromS 

Subscripts 

b  =  backward  reaction 
e  =  edge  condition 
/  =  forward  reaction 
w  =  wall  condition 

4.1  Introduction 

Stability  theory  of  compressible  boundary  layers  has  been  established  by  theoretical  studies 
by  Lees  and  Lin  [LL46],  Lees  and  Reshotko  [LR62],  and  Mack  [Mac65,  Mac69,  Mac75, 
Mac] .  More  recent  progress  in  the  understanding  of  stability  and  receptivity  of  hypersonic 
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boundary  layers  has  demonstrated  the  importance  of  instability  prehistory  analysis  [FKOl, 
FK02,  Fed03].  For  example,  hypersonic  instability  prehistory  reveals  the  following  features: 

1.  In  the  leading  edge  region,  discrete  modes  can  be  synchronized  with  fast  and  slow 
acoustic  modes.  This  synchronization  provides  an  efficient  channel  for  the  generation 
of  discrete  modes  by  acoustic  perturbations; 

2.  Further  downstream,  a  decaying  discrete  mode  can  be  synchronized  with  free-stream 
entropy  and  vorticity  perturbations.  This  indicates  that  entropy  and  vorticity  pertur¬ 
bations  can  generate  a  slightly  decaying  mode; 

3.  The  decaying  mode  also  can  be  synchronized  with  an  unstable  discrete  mode  (second 
Mack’s  mode)  leading  to  the  following  channel;  vorticity/entropy  decaying  discrete 
mode  second  mode. 

These  spectral  features  dehne  the  physics  of  boundary  layer  instabilities,  and  one  must 
explore  the  spectral  properties  in  detail  in  order  to  understand  the  receptivity  and  instability 
of  high-speed  boundary  layers. 

Almost  all  theoretical  and  computational  studies  of  the  stability  of  high-speed  boundary 
layers  have  used  the  calorically  perfect  gas  model.  Real  gas  effects  and  surface  chemistry 
bring  many  other  parameters  into  consideration  that  can  have  a  strong  impact  on  the  spec¬ 
trum  topology.  There  are  a  limited  number  of  publications  on  the  stability  of  boundary  layers 
including  real  gas  effects  [MA91,  Stu91,  SR94,  CVM97,  HCC97,  JSC98,  Ber98,  Lyt03,  LR05]. 
These  studies  provide  insight  into  real  gas  effects  in  boundary  layers  of  high-speed  air. 

Malik  and  Anderson  [MA91]  considered  chemical  equilibrium  boundary  layers  in  chemical 
equilibrium  high-speed  air  flow  past  a  flat  plate.  Real  gas  effects  were  introduced  through 
mean  flow  gas  properties,  while  the  stability  equations  remained  basically  the  same  as  in 
the  case  of  a  perfect  gas.  They  showed  that  real  gas  effects  lead  to  thinner  boundary  layers 
with  a  signihcantly  lower  maximum  temperature  than  in  the  calorically  perfect  gas  model. 
The  stability  analysis  indicated  that  real  gas  effects  stabilize  the  first  mode  and  destabilize 
the  second  one.  It  was  found  that  in  real  gases  the  peak  of  the  second  mode  is  shifted 
to  lower  frequencies.  Stuckert  [Stu91]  and  Stuckert  and  Reed  [SR94]  explored  the  effects 
of  nonequilibrium  chemical  reactions  in  air  on  the  stability  of  the  boundary  layer  over  a 
cone  at  high  Mach  numbers.  They  also  observed  the  frequency  shift  in  the  second  mode 
in  both  equilibrium  and  nonequilibrium  air  calculations.  Chang  et  al.[CVM97]  considered 
flows  in  chemical  equilibrium  and  nonequilibrium  using  the  parabolized  stability  equations. 
Hudson,  Chokani,  and  Candler  [HCC97]  extended  the  stability  analysis  to  flows  in  chemical 
and  thermal  nonequilibrium.  They  found  that  thermal  nonequilibrium  may  destabilize  the 
most  unstable  oblique  first  mode.  In  the  considered  examples,  they  did  not  observe  a  signif¬ 
icant  difference  in  the  second  mode  computations  using  models  of  chemical  equilibrium  or 
nonequilibrium  flows.  Johnson,  Seipp,  and  Candler  [JSC98]  explored  the  effect  of  the  free 
stream  total  enthalpy  on  the  stability  of  chemical  and  thermal  nonequilibrium  flows  of  air 
and  nitrogen.  They  showed  that  perturbations  of  concentrations  have  a  signihcant  impact  on 
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the  stability  of  high  enthalpy  air.  For  the  considered  cases,  they  suggested  that  the  total  en¬ 
thalpy  effect  in  air  is  due  to  oxygen  dissociation  and  that  the  dissociation  leads  to  absorption 
of  the  perturbation  energy.  Bertolotti  [Ber98]  studied  the  influence  of  rotational  and  vibra¬ 
tional  energy  relaxation  on  boundary-layer  stability.  He  showed  that  rotational  relaxation 
has  a  damping  effect  on  high-frequency  perturbations,  whereas  vibrational  relaxation  affects 
the  growth  rates  of  perturbations  by  changing  the  mean  flow.  Lyttle  [Lyt03]  and  Lyttle  and 
Reed  [LR05]  addressed  the  sensitivity  of  the  stability  analysis  to  the  thermodynamic  models, 
reaction  rates,  and  transport  models.  They  found  that  the  growth  rate  of  the  second  mode 
can  be  significantly  affected  by  the  chosen  model  of  the  chemically  reacting  air. 

In  order  to  assess  the  impact  of  real  gas  effects  on  stability  and  receptivity,  we  suggest 
building  simple  knowledge  blocks  based  on  the  analysis  of  such  gases  as  nitrogen  and  oxygen 
before  moving  on  to  the  analysis  of  air.  In  the  present  work  we  consider  the  inviscid  stability 
of  high-speed  boundary  layers  in  binary  mixtures  of  oxygen  and  nitrogen. 


4.2  The  real  gas  model 


We  consider  binary  mixtures  of  oxygen  and  nitrogen  assuming  that  they  are  in  local  thermal 
equilibrium;  that  is,  the  translational,  rotational  and  vibrational  degrees  of  freedom  are 
associated  with  a  single  temperature.  The  reaction  rates,  thermodynamic  properties  and 
transport  properties  are  modeled  following  [GYTL90]  .  The  mixtures  are  considered  to  be 
mixtures  of  perfect  gases.  The  enthalpy,  h,  and  frozen  specific  heat  at  constant  pressure, 
Cpf,  of  a  gas  mixture  can  be  determined  if  the  molar  fraction  of  each  species,  Xi,  is  known: 


h  ^  ^  Xi  hi , 
i 


i 


— 


P 


Here,  hi  and  Cp^  are  the  enthalpy  and  specific  heat  of  species  f,  respectively,and  {Ah{)T^^p 
is  the  heat  of  formation  of  species  i  at  temperature  T^e/.  Units  of  the  thermodynamic  and 
transport  properties  presented  in  this  section  are  specified  in  the  nomenclature.  [GYTL90] 
provides  curve  fits  for  hi  and  Cp^  as  functions  of  the  local  temperature  T. 

The  viscosity  of  the  gas  mixture  is  determined  by  the  approximate  formula 
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where  Mj  is  the  molecular  weight  of  species  i,  is  Avogadro’s  number,  and  Runiv  is  the 
universal  gas  constant.  The  average  collision  cross  section  hlT  ^  is  a  function  of  temperature 
(the  curve  £t  is  given  in  [GYTL90]). 

The  frozen  thermal  conductivity  Kf  of  the  mixture  is  calculated  according  to 


Kf  —  Rtr  +  Kint  (4-3) 

where  Kt^.  is  the  translational  component  of  the  thermal  conductivity  and  Kint  is  the  compo¬ 
nent  of  thermal  conductivity  associated  with  the  internal  excitation  energy  of  the  molecules. 
They  are  given  by 
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where  k  is  the  Boltzmann  constant.  Quantities  Gg’^^  and  Ag^  are  functions  of  temperature. 
Their  curve  fits  are  given  in  [GYTL90]. 

The  binary  diffusion  coefficient  Dij  is  evaluated  according  to  the  following  equation: 
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where  Ag^  is  defined  in  Eq.  (4.4). 

The  rate  of  oxygen  dissociation  is  independent  of  the  third  body  M^at  participating  in 
the  collision.  The  oxygen  dissociation  reaction  O2  +  Meat  ^  20  -|-  Meat  is  associated  with  the 
following  forward  and  backward  reaction  rates  k!"p  and  A;g\  respectively: 


k^P  =  3.61  X  10^®T-^-°exp  (-5.94  x  lOV^)  ,  (4.6) 

A;g^  =  3.01  X  10^^T-°-^ 

The  nitrogen  dissociation  rate  depends  on  the  third  body.  The  reaction  A2  -|-  X2  ^  2N  +  N2 
has  the  reaction  rates 
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(4.7) 


kf  =  1.92  X  10^^T-°-®exp  (-1.131  x  lOV^)  , 

=  1.09  X  10^®T-°■^ 

while  the  reaction  iV2  +  iV  ^  2N  +  N  has  the  reaction  rates 

kf  =  4.15  X  lO^^T-^-^exp  (-1.131  x  lOV^)  , 

=  2.32  X  X  K,or,  (4.8) 

Kcor  =  1.0155. 

We  introduce  the  correction  factor,  K^or  ,  in  Eq-  (4.8)  in  order  to  have  the  same  equilibrium 
constant  for  nitrogen  dissociation  when  the  third  body  is  N2  or  N,  whereas  the  reaction 
rates  for  nitrogen  in  [GYTL90]  give  different  ratios  kj  /k^  and  kj  /k^  . 


4.3  The  mean  flow  problem 


4.3.1  Governing  equations 

The  two-dimensional  mean  flow  profiles  have  been  found  numerically  by  solving  the  boundary- 
layer  equations  for  a  reacting  gas  mixture  [Dor62]: 
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Here,  Pe  is  the  pressure  at  the  edge  of  the  boundary  layer;  /  =  h  +  u^ / 2  is  the  total  enthalpy; 
Pr  and  Le  are  the  Prandtl  and  Lewis  numbers,  respectively;  ZI12  is  the  binary  diffusion 
coefficient;  q  is  the  mass  fraction  of  species  i;  and  Wi  is  the  net  mass  rate  of  production  of 
species  i.  The  equation  of  state  for  the  gas  mixture  is 


P  (4.10) 

i 

where  Ri  =  Runiv/Mi  is  the  gas  constant  of  species  i. 

In  the  case  of  a  binary  mixture  of  oxygen,  the  net  mass  production  rate  Wi  in  Eq.  (4.10) 
is  given  by 
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W,  =  2k'/yM,^  +  2k<‘yM2-^  (^^  +  , 

In  Eq.  (4.11)  and  what  follows,  we  use  subscripts  1  and  2  for  molecules  and  atoms,  respec¬ 
tively.  For  a  binary  mixture  of  nitrogen  we  have 
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Rather  than  using  the  continuity  equations  for  both  species,  we  apply  the  total  mass 
conservation  law  and  use  only  the  continuity  equation  for  atoms  (or  molecules). 

The  equations  are  transformed  using  the  Levy-Lees-Dorodnitsyn  transformation  [Dor62]: 
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0  0 

where  subscript  e  indicates  the  value  at  the  edge  of  the  boundary  layer.  Using  the  stream 
function  ip  {P,,p),  one  can  End  the  streamwise  velocity  component 


u  =  Uef  (C,h) , 

where  ’prime’  indicates  the  partial  derivative  with  respect  to  77. 
The  transformed  boundary  layer  equations  are  written  as 


(4,14) 


{Cf")'  +  ff"  +  PH  'p-if  f  ^  (f'%;  -  4%)  • 

C*  ^  f  C  f  1  ^  ^  hiCie  f  ul\  f  1  T  ^  ^ 

7;  - 1 )  1.  — 9.  +1-  (  pc-qc'// 


Sm  \  Le 


-44-44 


le  \\Pr 


_  fyL  2^d^ 

V(9^  \Sm^y  Cie  d^  ppeUjpeCy 


(4.15) 


61 


2^dUe  ,  .  .  .  Pll 

PH  =  — 9  [i,  V)  =  I/h:  Qi  =  Ci/Cie,  C  =  - , 

tie  Pef^e 

where  Sm  =  p/pDi2  is  the  Schmidt  number.  Equations  (4.15)  are  to  be  solved  subject  the 
following  boundary  conditions: 


h  =  0  :/  =  /'  =  0;  (4.16) 

p  ^  oo  :  f  ^  l,g  ^  l,qi  ^  1. 

In  addition  to  (4.16),  we  have  to  formulate  boundary  conditions  at  the  wall  for  the  species 
mass  fraction,  q,  and  the  enthalpy  (or  temperature),  g.  The  boundary  condition  on  the 
species  mass  fraction  depends  on  the  recombination  rate  at  the  wall  [Dor62] .  The  boundary 
condition  on  the  enthalpy  depends  on  the  specihc  case  under  consideration.  It  can  be  de¬ 
termined  by  specifying  the  wall  temperature  or  the  heat  flux.  In  the  case  of  the  wall  heat 
flux  being  specihed,  one  must  take  into  account  that  the  heat  flux  has  contributions  from 
thermal  conductivity  and  the  species  diffusion  toward  the  wall  [GYTL90,  Dor62]. 

4.3.2  Numerical  method 

The  computational  domain  in  ^  and  p  is  discretized,  and  a  second-order  hnite-difference 
approximation  for  Eqs.  (4.15)  similar  to  that  found  in  [Ceb99]  is  utilized.  The  grid  points 
are  dehned  as 


=  0,  r  =  r-'  +  kn,n  =  1,  2, ...,  N;  (4.17) 

Vo  =  0,  Pj  =  Pj-i  +  hj,j  =  1,  2, ...,  J, 

where  hj  stands  for  the  grid  step  size  (not  to  be  confused  with  enthalpy  in  sections  2  and  3). 
For  the  binary  mixture,  we  use 


CleQl  +  C2eQ'2  —  1 

Introducing  the  new  dependent  variables 

qi=Qi,  q[  =Q2j  =  Fi,  f  =  F2,  r  =  Fs,g  =  G,,g'  =  G2,  (4.18) 

the  first-order  equations  are  discretized  and  written  for  the  midpoint  (^^'^,pj-i/2)  ■ 
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Equations  (4.15)  are  also  discretized  at  the  midpoint  77j_i/2)  .  The  numerical  algo¬ 

rithm  is  organized  by  solving  the  momentum  equation  and  the  equation  of  species  conserva¬ 
tion  separately  from  the  energy  equation.  That  is,  initial  values  for  Gi  and  G2  are  chosen 
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and  the  five  equations  representing  momentum  and  mass  conservation  are  solved  to  obtain 
Fi,  F2,  F3,  Qi,  and  Q2  based  on  the  given  Gi  and  G2  ■  New  values  Gi  and  G2  are  than 
found  by  solving  the  two  equations  representing  energy  conservation  based  on  the  values  for 
Fi,  F2,  F3,  Qi,  and  Q2  generated  in  the  previous  step.  The  new  Gi  and  G2  are  than  used  to 
obtain  updated  values  of  Fi,  F2,  F3,  Qi,  and  Q2  from  the  momentum  and  mass  conservation 
equations.  This  iterative  process  is  continued  until  a  solution  within  a  certain  tolerance  is 
reached. 

The  nonlinear  algebraic  system  of  equations  is  solved  using  Newton’s  iterations  that  lead 
to  a  linear  system  of  equations  for  the  ’’perturbation”  of  each  function  : 

=  Qh-b +  (4.20) 

where  index  i  indicates  the  iteration  number.  The  algorithm  leads  to  a  block  tridiagonal 
structure  of  equations  that  can  be  expressed  in  matrix-vector  form 


A5  =  f,  (4.21) 

where  the  coefficient  matrix  A  and  the  vector  r'are  known.  For  example,  Newton’s  iterations 
for  the  discretized  equations  for  SFi,SF2,SF3,SQi,  and  <5^2  lead  to  the  following  structure 
of  the  algebraic  system: 
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where  the  vector  r  is  known  from  the  previous  iteration.  MatricesA^,  Fj,  and  Gj  are  5x5 
matrices.  Their  explicit  form  is  given  in  [UlklO].  The  linear  set  of  equations  (4.21)  is  solved 
using  the  block  elimination  method  [Ceb99]. 


4.3.3  Some  examples  of  the  mean  flow  profiles 

For  the  purpose  of  testing  the  boundary  layer  solver,  we  have  considered  examples  corre¬ 
sponding  to  the  first  and  second  Probstein’s  integrals  [Dor62]. 

If  die/dx  =  0,  dcie/d^  =  0  ,  Pr  =  Sm  =  1  ,  and  W  =  0  (frozen  flow),  the  first  Probstein’s 
integral  is 
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Flat  Plate.  0+0,.  U  =2  km/s,  P  =0.01  atm,  T  =3000  K 

2  e  e  e 


Figure  4.1:  Comparison  of  numerical  results  for  mass  fraction  with  the  hrst  Probstein  inte¬ 
gral. 


Qi  =  (9  -  9w) ,  (4.23) 

-L  9w 

where  subscript  w  stands  for  values  at  the  wall.  Figure  4.1  shows  a  comparison  of  the 
numerical  results  with  the  solution  given  by  Eq.  (4.23)  for  a  binary  mixture  of  oxygen 
in  the  boundary  layer  past  a  flat  plate  at  the  edge  velocity,  pressure  and  temperature  of 
Ue  =  2km/s  ,  Pe  =  O.Olatm  ,  and  Te  =  3000K,  respectively.  The  edge  mass  fractions  were 
chosen  according  to  the  equilibrium  state  of  the  mixture,  and  =  0.8.  The  length  scale  is 
H  =  V^/{PeUe). 

If  due/dx  =  0  ,  dcie/d^  =  0  ,  Sm  =  1  ,  and  W  =  0  (frozen  flow),  the  second  Probstein’s 
integral  is 

Qi  =  (?i)^  +  [1  -  (?i)^]  /'•  (4.24) 

A  comparison  of  the  numerical  results  with  the  second  Probstein’s  integral  is  shown  in  Fig. 
4.2.  One  can  see  good  agreement  between  the  numerical  results  and  Probstein’s  integrals. 

Another  test  is  illustrated  by  Figs.  4.3  and  4.4,  where  the  self-similar  solution  was 
computed  by  solving  a  system  of  ordinary  differential  equations  (4.20)  for  a  boundary  layer 
past  a  flat  plate  at  Mach  number  M  =  2  ,  edge  temperature  =  300K  ,  and  Prandtl  number 
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Flat  Plate.  0+0,.  U  =2  km/s,  P  =0.01  atm,  T  =3000  K 

2  e  e  e 


Figure  4.2:  Comparison  of  numerical  results  for  mass  fraction  and  the  second  Probstein’s 
integral. 
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Flat  Plate.  0+02.  P  =1  atm,  T  =  300  K,  U  =  660.7  m/s  (M=2) 

e  e  e  '  ' 


Figure  4.3:  Comparison  of  the  velocity  profile  with  the  self-similar  solution. 


Pr  =  0.7.  In  order  to  compare  the  results  from  the  new  solver  with  the  self-similar  solution  for 
a  non-reacting  calorically  perfect  gas,  we  adjusted  the  formulation  in  the  new  solver  to  have 
a  constant  Prandtl  number  and  to  consider  a  gas  with  constant  specific  heat.  Additionally, 
we  let  the  viscosity  be  given  by  Sutherland’s  viscosity  as  a  function  of  temperature. 


4.4  Inviscid  stability  analysis.  Governing  equations 
and  numerical  method 


4.4.1  Governing  equations 

We  consider  three-dimensional  perturbations  in  a  two-dimensional  boundary  layer  flow  using 
the  quasi-parallel  flow  approximation.  The  linearized  governing  equations  in  the  inviscid 
limit  can  be  written  as  follows: 


du'  du'  ,du  dp' 


dv'  dv'  dp' 

P^+  =  0, 

dt  dx  dy 


(4.25) 

(4.26) 
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Figure  4.4: 


Flat  Plate.  0+02.  P  =1  atm,  T  =  300  K,  U  =  660.7  m/s  (M=2) 

e  e  e  '  ' 


Comparison  of  the  temperature  profile  with  the  self-similar  solution. 
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dw'  dw'  dp' 

fdt  +'’"ax  =°' 

(4.27) 

dp'  dp'  du'  dv'  fdp  dw' 

m  ^f’dx  ^f’dy  dy^f’  dz 

(4.28) 

dV  dp'  dl'  .  dl 

(4.29) 

'  dc' 

^  /•  1  o\ 

'’at  = 

(4.30) 

where  ‘prime’  indicates  a  perturbation.  For  a  binary  mixture,  c[  +  C2  =  0  ,  one  can  find 


,,,,  ,  ,dm  \  ,  ^(dw,  ,  dW,  ]  , 

=  PCi  I  wn - ET-  1  +  P  I  ir-ci  +  yr-C2  1  +  , 


dpi  dp2  J  y  dpi  dp2 

p  =  c[pT  {Ri  -  R2)  +  p'RT  +  pRT' , 


dT 


(4.31) 

(4.32) 


I'  —  c' I  {hi  —  h2)  +  T'Cpf  +  uu\ 
fdhi  dh2\ 

Cp  f  —  Cl  ——  +  C2  ^777  . 


(4.33) 


dT  ^  dT 


The  solution  of  the  linearized  equations  is  considered  in  the  form  of  normal  modes 


q'  (x,  y,  z,t)  =  q  {y)  exp  [i  (ax  +  fdz  -  ut)] , 
and  we  arrive  at  the  following  system  of  ordinary  differential  equations: 

/  -  /I  -  \  ^dtp 

I  (au  —  oj)  p  +  ip  (au  +  jjw)  +  p- — h  v—  =  0, 

dy  dy 

^du 

ip  (au  —  uj)u  +  pv- — h  lap  =  0, 
dy 

s  ^  dp 

ip  (au  —  a;)  n  +  —  =  0, 
dy 

ip  (au  —  uj)w  +  i(3p  =  0, 

\  f  - 

ip  (au  —  ui)  I  +  lup  +  pv—  =  0, 

dy 

ip  (au  —  a;)  Cl  +  pv-r-  =  ciAi  +  pBi  +  CiT, 
dy 


(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 

(4.39) 

(4.40) 
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Ai  =  p 


dWi  dWi\  _(dWi  dWi  \  _  dWi 

Similar  to  the  inviscid  stability  analysis  in  non-reacting  gases,  one  can  derive  from  Eqs. 
(4.35)  -  (4.40)  a  single  eqnation  for  the  pressnre  pertnrbation: 


ay  dy^ 


(4.41) 


n  =  -  (a^ 


+  /3") 


(a.u—uj)  I  {Ri—R2)G\  (olu—uj) 

T  ^  R  Wi 


(i?i  —  R‘2)  pBi  {au  —  uY 
R  ~Di 

[Di  -  (hi  -  hY)] 


Fo  = 


[{hi  —  ^2)  Gi  +  DiCpf] 

2a  du  (i?i  —  R2)  dci 

1 - V 1 - 7T~5 — 

[au  —  00)  dy  DiR  dy 


1  Gii?i 
T  “  DiR 


f  uf\D,-ip(au-uj)(h,-h,)1^ 


,  du 


dy  dy 


dci 


dy 


[(hi  —  h2)  Gi  +  DiGpf] 


Di  =  ip  {au  -  cn)  p^^^  Bi  -  Ai, 

it 

G,  =C, -S,^, 

where  Ai,Bi,  and  Gi  are  dehned  by  Eq.  (4.40). 

Also  similar  to  the  stability  analysis  of  non-reacting  gases,  for  the  nnmerical  solving  of 
the  stability  eqnations  we  nse  a  system  of  two  differential  eqnations  written  in  dimensionless 
form: 


^  =  Fuv,  (4.42) 

dy 

5  =  F21P  +  F22V, 
dy 

where  the  pressnre  pertnrbation  is  scaled  by  peul  ,  and  the  velocity  pertnrbation  is  scaled 
by  Me  .  The  dimensionless  coordinate  y  =  y/5  is  dehned  with  the  help  of  a  length  scale  S 
(nsnally  we  nse  the  length  scale  S  =  H  =  y^/  {Pe^e)  or  h  =  Flu  =  v^/  (Pe'^e))-  Coefficients 
Fi2  ,  ^"21  ,  and  F22  are  given  by 


T"i2  — 


Pe 
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r2^  =  -I  \  ao - 


Ue 


Up 


PeUl  {Ri-R2)^  Te^  Pe  (M)'  +  W) 


Pe 


„  -Qi  -  - 

R  T  p 


a6^  -  ^ 

Re  Re 


F„  =  *  I  ai  “  - 

UeJ 


Up 


{Ri 


iaS 


du/up 


+  i~ 


pe 

P 


dp/ Pe 


_  ^11022  ~  012^21  ^  _ 

Vl  —  - )V2  — 

011022  ~  O12O21 

_  &21O11  —  O21&II  „  _ 

Oi  —  - ,02  — 

O11O22  ~  O12O21 


ad^-^Y  ^y 

Re  Re  J  \  J 

&12O22  ~  O12&22 


u  uj6\  6Bi{Ri—R2 

On  =  M  ad - + - ^ - - - 

Up  Up  J  Up  R 


011O22  ~  012O21 
&22O11  —  021&12 

5 

CillCi22  ~  (^12(^21 

Pe  dAi 


012  — 


PedTpC^  ,  5B,R 


P  Pe'^e 


U 

Up 


+  -  (  ad- —  —  )  (  -r  ~  ~r 


P  ^epe 

ujd\  f  hr 
Up  I  \  L  L 


u  cc;(5\  TeC  f 
022  =  ad - — 


611  = 


Up  Ue  J 

Peul  dB 


Pe 


1  dci 

-M2  = 

dy 


U  ud\  PeU 


621  =  ad - 


Up 


Up 


'-M2  =  i 


.fd{I/Ie)  ul  u  d^^ 


pie  V  dy  Ip  Ue  dy 

In  order  to  arrive  at  the  eigenvalue  problem,  one  must  formulate  homogeneous  boundary 
conditions  for  Eqs.  (4.42) 


p  — 00  :  p  — 0;  (4.43) 

p  =  0  :  h  =  0  (or  dp/dy  =  0). 

4.4.2  Numerical  method 

In  the  present  work  we  consider  the  stability  of  boundary  layers  in  the  temporal  framework. 
In  this  case,  a  and  (I  are  real  parameters  and  the  complex  frequency  uj  is  found  as  an 
eigenvalue.  Inviscid  analysis  requires  using  a  path  of  integration  in  the  complex  plane  y  that 
is  below  the  critical  point  pc  defined  by  the  equation  of  the  critical  layer  [Mac65]: 

u  (ye)  =  c. 
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where  c  =  loja  is  the  complex  phase  velocity.  In  the  numerical  stability  analysis  of  the 
inviscid  equations  for  perturbations  in  non-reacting  boundary  layer  flow,  Mack  [Mac65]  used 
four  terms  in  the  Taylor  series  in  order  to  determine  the  velocity  and  temperature  prohles 
of  the  mean  flow  at  complex  y^.  The  same  approach  was  used  in  [FTIO].  For  the  present 
analysis  of  reacting  real  gas  flow,  this  approach  is  not  realistic  because  of  the  complexity  of 
the  equations  for  the  mean  flow.  In  addition,  many  thermodynamic  properties  of  the  mixture 
are  determined  using  curve  htting,  which  prevents  hnding  higher  derivatives  with  acceptable 
accuracy.  Therefore,  we  restrict  ourselves  by  considering  only  unstable  perturbations  (c*  = 
Im(c)  >  0)  and  using  the  integration  path  along  the  real  axis  y. 

We  use  the  fourth-order  Runge-Kutta  method  with  a  constant  step  size  to  solve  Eqs. 
(4.42)  starting  outside  the  boundary  layer  at  ymax  and  integrating  toward  the  wall.  Outside 
the  boundary  layer,  the  coefficient  F22  in  Eqs.  (4.42)  is  equal  to  zero,  and  one  can  End 

V  (j/)  ~  exp  (Aj/) , 

A  =  \/ F\2  iVniax)  E21  (l/  max  ),Re(A)<0. 

For  the  purpose  of  code  testing,  we  reproduce  results  from  [FTIO]  for  a  non-reacting, 
calorically  perfect  gas.  Figure  4.5  illustrates  a  comparison  for  a  two-dimensional  perturbation 
(/5  =  0)  .  In  order  to  compare  with  the  results  of  [FTIO],  we  must  enforce  properties  of  a 
calorically  perfect  gas  on  the  coefficients  of  Eqs.  (4.42). 

4.5  Inviscid  stability  analysis.  Numerical  examples  for 
reacting  binary  mixtures 

4.5.1  Oxygen 

In  the  present  work,  we  are  interested  in  studying  the  impact  of  real  gas  effects  on  flows  in 
thermal  equilibrium.  In  this  sense,  the  work  is  aligned  with  [MA91].  Rather  than  solving 
a  complete  set  of  ordinary  differential  equations  including  thermal  conductivity,  diffusivity, 
and  viscosity  for  perturbations  in  a  boundary  layer  of  equilibrium  air,  we  begin  with  an 
inviscid  analysis  and  consider  binary  mixtures  of  chemically  reacting  boundary- layer  flows. 
In  all  examples,  we  consider  two-dimensional  perturbations  with  (3  =  0. 

The  first  example  is  the  equilibrium  free  stream  flow  of  oxygen  at  zero  pressure  gradient 
at  pressure  Pe  =  0.01  atm,  edge  velocity  Ue  =  4.273  km/s,  and  edge  temperature  Tg  =  2000 
K.  The  edge  parameters  in  the  binary  mixture  of  oxygen  correspond  to  the  frozen  Mach 
number  M  =  5.2.  In  this  example,  we  consider  an  adiabatic  and  non-catalytic  wall.  The 
plate  length  is  equal  to  2  m,  and  the  total  enthalpy  of  the  free-stream  flow  is  11  MJ/kg.  The 
total  enthalpy  is  approximately  the  same  as  in  [JSC98]  (9.31  MJ/kg)  where  the  boundary- 
layer  flow  of  air  was  considered  and  a  total  enthalpy  effect  was  observed.  Figures  4.6  and 
4.7  illustrate  a  comparison  of  the  velocity  and  temperature  profiles  of  the  reacting  gas  with 
the  self-similar  solution  for  non-reacting  gas  flow  with  specific  heat  ratio  1.4.  In  agreement 
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Flat  Plate.  Adiabatic  wall.  M  =  3.8 


Figure  4.5:  Comparison  of  the  growth  rates  of  the  present  work  and  Ref.  [FTIO]  for  a 
non-reacting  calorically  perfect  gas. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.6:  Comparison  of  the  velocity  profiles  in  real  and  perfect  gases. 


with  observations  of  [MA91],  one  can  see  that  the  thickness  of  the  boundary  layer  in  a  real 
gas  is  smaller  and  that  there  is  a  dramatic  effect  on  the  temperature  prohle. 

In  order  to  examine  the  influence  of  real  gas  effects  on  the  stability  equations,  we  consider 
coefficients  Fi  and  F2  in  Eq.  (4.41).  In  the  absence  of  a  mass  fraction  perturbation  and 
assuming  a  calorically  perfect  gas,  Fi  has  the  following  limit: 

J-  rjn  * 

where  U  =  u/ue  and  T  =  T /T^.  Figure  4.8  shows  the  real  part  of  the  coefficient  Fi  for  a 
reacting  gas,  for  a  calorically  perfect  gas  using  Eq.  (4.44)  with  the  velocity  and  temperature 
prohles  for  the  real  gas  mean  flow,  and  according  to  Eq.  (4.44)  with  the  velocity  and 
temperature  prohles  for  a  non-reacting  gas.  Although  there  are  some  differences  depending 
on  the  gas  model  used  for  the  coefficient  Fi,  the  curves  have  the  same  qualitative  behavior. 
To  evaluate  the  impact  of  real  gas  effects  on  the  coefficient  F2  in  Eq.  (4.41),  we  consider 


2  dU 


(4.45) 


In  the  case  of  a  non- reacting  and  calorically  perfect  gas,  F^  =  —T'/T.  Figure  4.9  illustrates 
the  difference  between  reacting  and  non-reacting  gases  for  the  coefficient  F^.  Qualitatively, 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.7:  Comparison  of  nthe  temperature  profiles  in  real  and  perfect  gases. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.8:  Effect  of  the  real  gas  flow  on  the  coefficient  Fi  on  Eq.  (4.41). 


the  behavior  of  the  coefficient  is  the  same  for  both  gas  models.  These  comparisons  indicate 
that  the  structure  of  the  spectrum  should  be  qualitatively  the  same. 

Figure  4.10  shows  the  imaginary  part  Cj  versus  a.  One  case  shown  in  Fig.  10  corresponds 
to  integration  with  4000  intervals  in  the  boundary  layer  on  the  interval  [0,14.23];  the  other 
case  corresponds  to  8000  intervals.  This  illustrates  that  the  result  is  independent  of  the  step 
size  used  in  integration.  One  can  see  the  typical  structure  of  the  spectrum  having  the  hrst 
and  the  second  modes  (the  reader  can  End  details  regarding  the  terminology  in  [FTIO]). 

Figure  4.11  shows  a  comparison  of  c*  for  a  real  gas  and  a  calorically  perfect  gas.  One 
can  see  that  real  gas  effects  stabilize  the  hrst  mode  and  destabilize  the  second  mode.  This 
phenomenon  was  noted  in  [MA91].  It  was  interpreted  as  the  result  of  wall  cooling  due  to 
real  gas  effects.  However,  interpretation  of  the  result  as  a  pure  temperature  effect  might  be 
incorrect  as  we  will  see  in  our  discussion  of  the  results  for  nitrogen  in  the  following  section. 

The  typical  eigenfunctions  for  mode  1  and  mode  2  are  shown  in  Figs. 4. 12  and  4.13  as 
the  magnitude  of  the  pressure  perturbation.  The  shapes  of  these  eigenfunctions  are  similar 
to  those  observed  for  non-reacting  gases  [Mac69] . 

The  real  gas  effects  appear  in  the  governing  stability  equations  as  non-calorically  perfect 
gas  properties,  as  changes  in  the  mean  velocity  and  temperature  profiles,  and  as  the  presence 
of  perturbations  of  the  species  mass  fractions.  If  we  exclude  the  equation  for  mass  fraction 
perturbations  from  the  governing  equations,  the  eigenvalues  are  changed  as  shown  in  Fig. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.9:  Effect  of  the  real  gas  flow  on  the  coefficient  on  Eq.  (4.45). 


76 


dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.10:  Comparison  of  Cj  for  two  different  step  sizes  used  in  integration. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.11:  Comparison  of  c*  for  real  gas  and  calorically  perfect  gas. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 
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Figure  4.12:  Pressure  distribution  across  the  boundary  layer  for  mode  1. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 
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Figure  4.13:  Pressure  distribution  across  the  boundary  layer  for  mode  2. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.14:  Pressure  distribution  across  the  boundary  layer  for  mode  2. 


4.14.  However,  the  spectrum  structure  remains  qualitatively  the  same.  This  result  agrees 
with  studies  of  the  free-stream  total  enthalpy  effect  in  [JSC98].  In  addition  to  modes  1 
and  2  that  are  subsonic  (the  external  flow  in  the  frame  moving  with  the  phase  speed  of 
the  perturbations  is  subsonic),  Fig.  4.14  also  shows  two  supersonic  perturbations  found  in 
this  study.  As  was  pointed  out  by  Mack  [Mac69] ,  the  supersonic  perturbations  have  a  much 
smaller  growth  rate  than  the  second  mode  has.  Chang  et  ah  [CVM97]  also  found  unstable 
supersonic  discrete  modes  in  boundary-layer  flows  of  air  within  the  parabolized  stability 
equations  framework. 

The  effect  of  mass  fraction  perturbations  on  the  stability  of  the  binary  mixture  of  oxygen 
in  Fig.  4.14  was  obtained  within  the  temporal  stability  framework.  In  the  spatial  framework, 
we  observe  the  same  phenomenon  (see  Fig.  4.15  where  the  imaginary  part  of  the  streamwise 
wave  numbers  is  shown  as  a  function  of  frequency). 

In  the  next  example,  we  consider  the  effect  of  wall  cooling  on  the  inviscid  two-dimensional 
perturbations.  The  wall  temperature  is  =  1000  K,  and  the  mass  fraction  of  atoms  on  the 
wall  is  C2w  =  0.  One  can  see  from  Fig.  4.16  that  cooling  destabilizes  the  second  mode  as  is 
known  from  studies  of  non-reacting  gases.  There  is  not  a  strong  effect  on  the  hrst  mode  at 
the  chosen  flow  parameters. 
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dPe/dx=0.  Oxygen.  Adiabatic+non-catalitic  wall 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 
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Figure  4.15:  Pressure  distribution  across  the  boundary  layer  for  mode  2. 
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dPe/dx=0.  Oxygen. 

Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s  (M^  =  5).  Equilibrium  free  stream. 


Figure  4.16:  Pressure  distribution  across  the  boundary  layer  for  mode  2. 
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dPe/dx=0.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  11  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


Figure  4.17:  Temperature  profiles  in  binary  mixtures  of  oxygen  and  nitrogen. 

4.5.2  Nitrogen 

In  order  to  compare  stability  characteristics  of  nitrogen  and  oxygen,  we  consider  an  equi¬ 
librium  free  stream  flow  of  a  binary  mixture  of  nitrogen  at  the  same  conditions  chosen  for 
oxygen:  pressure  Pe  =  0.01  atm,  edge  velocity  Ue  =  4.273  km/s,  and  edge  temperature 
Tg  =  2000  K.  The  edge  parameters  correspond  to  the  frozen  Mach  number  M  =  4.9.  In 
order  to  explain  the  differences  in  these  two  examples  (binary  mixture  of  oxygen  and  bi¬ 
nary  mixture  of  nitrogen).  Figs.  4.17  and  4.18  show  the  temperature  and  mass  fraction  of 
molecules  in  the  boundary  layers.  One  can  see  that  the  temperature  in  the  oxygen  mixture 
is  lower  than  in  the  binary  mixture  of  nitrogen  due  to  dissociation. 

One  could  expect  that  the  lower  temperature  of  the  wall  should  destabilize  the  second 
mode  as  was  illustrated  in  Fig.  4.11  by  comparison  of  the  results  for  a  perfect  gas  with 
results  for  a  binary  mixture.  However,  in  the  present  case,  the  second  mode  in  nitrogen  is 
more  unstable  (see  Fig.  4.19).  This  means  that  an  interpretation  of  the  real  gas  effects  on 
the  second  mode  may  be  more  complicated  than  simply  a  pure  temperature  effect.  In  the 
present  example,  prohles  of  specihc  heats  of  oxygen  and  nitrogen  in  the  boundary  layers  are 
qualitatively  different.  Figure  4.20  shows  the  distribution  of  the  frozen  specihc  heats,  Cp/, 
in  boundary  layers  of  oxygen  and  nitrogen.  One  can  see  the  qualitative  difference:  while  the 
specihc  heat  of  nitrogen  increases  in  the  vicinity  of  the  wall  due  to  the  temperature  ehect. 


84 


Mass  Fraction  (Molecules) 


dPe/dx=0.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =11  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


y/H 


Figure  4.18:  Molecular  mass  fractions  in  binary  mixtures  of  oxygen  and  nitrogen. 
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dPe/dx=0.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  11  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


Figure  4.19:  Comparison  of  Cj  for  boundary  layers  in  nitrogen  and  oxygen. 


the  specific  heat  of  oxygen  is  affected  by  dissociation  near  the  wall.  This  could  be  the  reason 
for  the  stabilizing  effect  on  the  binary  mixture  of  oxygen  in  comparison  with  the  binary 
mixture  of  nitrogen. 

Because  there  is  no  signihcant  dissociation  in  nitrogen,  there  is  no  difference  between  the 
results  with  and  without  mass  fraction  perturbations  in  the  stability  equations  (Fig.  4.21). 

Figures  4.22  and  4.23  illustrate  coefficients  Fi  and  F3  in  Eqs.  (4.41)  and  (4.45),  respec¬ 
tively,  for  the  case  of  nitrogen  at  a  =  0.23,  c  =  0.914  -|-  ffi.023.  One  can  see  that  they  have 
the  same  qualitative  behavior  as  in  the  case  of  oxygen  presented  in  Figs.  4.8  and  4.9. 


4.6  Conclusion 

The  considered  examples  of  boundary  layer  flows  in  binary  mixtures  of  oxygen  and  nitrogen 
indicate  that  the  real  gas  effects  stabilize  the  hrst  mode  and  destabilize  the  second  one  in  the 
inviscid  limit  of  the  stability  equations  (see  [FTIO]  regarding  the  terminology).  The  result 
is  in  agreement  with  [MA91]  where  this  observation  was  made  for  boundary  layer  flow  for 
air  in  chemical  and  thermal  equilibrium  with  the  stability  equations  including  the  effect  of 
viscosity,  diffusivity,  and  thermal  conductivity.  Real  gas  effects  lead  to  a  decrease  in  wall 
temperature  and  can  have  an  impact  on  the  perturbations  similar  to  the  wall  cooling  effect 
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dPe/dx=0.  Nitrogen.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  1 1  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


y/H, 


Figure  4.20:  Comparison  of  frozen  specific  heats  in  the  binary  mixtures  of  oxygen  and 
nitrogen. 
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dPe/dx=0.  Nitrogen.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  1 1  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


Figure  4.21:  Comparison  of  c*  for  boundary  layer  in  nitrogen  when  perturbations  of 
fractions  are  included/excluded  in  the  stability  equations. 


dPe/dx=0.  Nitrogen.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  1 1  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


y/H 

Figure  4.22:  Coefficient  Fi  in  Eq.  (4.41).  Nitrogen,  a  =  0.23,  q  =  0.914  +  i0.023. 
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Re[F3] 


dPe/dx=0.  Nitrogen.  Adiabatic+non-catalitic  wall.  Total  Enthalpy  =  1 1  MJ/kg 
Pe=0.01  atm,  Te  =  2000  K,  Ue=4.273  km/s.  Equilibrium  free  stream. 


y/H 

Figure  4.23:  Coefficient  in  Eq.  (4.45).  Nitrogen,  a  =  0.23,  q  =  0.914  +  i0.023. 
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in  boundary  layers  of  non-reacting,  calorically  perfect  gases.  However,  such  an  interpretation 
of  the  wall  temperature  effect  might  be  misleading  when  binary  mixtures  of  nitrogen  and 
oxygen  are  compared.  In  the  considered  example,  nitrogen  does  not  experience  signihcant 
dissociation,  whereas  oxygen  goes  through  dissociation.  The  difference  leads  to  qualitative 
differences  in  the  specihc  heats,  and  the  second  mode  in  nitrogen  is  more  unstable  than  the 
second  mode  in  oxygen  in  spite  of  the  higher  wall  temperature. 

The  analysis  of  the  coefficients  in  the  Eq.  (4.41)  for  pressure  perturbations  in  the  inviscid 
limit  indicates  that  the  coefficients  qualitatively  have  the  same  structure  as  in  the  case  of  non¬ 
reacting  calorically  perfect  gases.  Thus,  the  conclusions  about  multiple  unstable  eigenmodes 
in  high-speed  boundary  layers  should  be  applicable  to  the  considered  examples  of  reacting 
boundary  layers. 

The  inviscid  stability  analysis  can  be  helpful  in  understanding  the  structure  of  the  spec¬ 
trum  due  to  the  simplicity  of  the  governing  equations.  Because  there  is  a  limited  number  of 
studies  examining  reacting  boundary  layers,  it  should  be  useful  to  initially  explore  the  real 
gas  effects  in  binary  mixtures  of  oxygen  and  nitrogen  before  moving  on  to  the  hve  species 
model  of  air.  The  present  work  demonstrates  only  a  small  sampling  of  examples  of  the  flow 
conditions,  and  studies  considering  a  broader  range  of  flow  parameters  are  still  required. 
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Chapter  5 

Flow  instabilities  and  transition 


This  chapter  is  based  on  the  article  [TumlO] 

5.1  Introduction 

The  hrst  study  of  laminar-turbulent  transition  was  carried  out  by  [Rey83]  who  attributed 
the  phenomenon  to  flow  instability.  The  importance  of  the  problem  for  design  of  laminar 
wings  for  commercial  aircraft,  prediction  of  viscous  drag  and  heating  of  the  surface  of  a 
hypersonic  vehicle  and  for  many  other  technical  applications  stimulated  broad  theoretical 
and  experimental  studies  of  flow  instabilities.  Today,  one  can  hud  a  vast  bibliography  on 
stability  and  transition  of  shear  flows  ([DR81,  SHOl,  CJJ67]).  There  are  a  variety  of  paths 
to  turbulence  depending  on  the  level  of  free-stream  disturbances  and  also  the  nature  of 
disturbances  introduced  into  the  flow  from  a  wall  (for  example,  isolated  and  distributed 
roughness).  There  are  also  different  mechanisms  governing  the  flow  instability.  In  order  to 
develop  and  use  a  reliable  tool  for  transition  prediction,  one  has  to  be  familiar  with  the  broad 
picture  of  the  transition  phenomenon.  A  short  article  in  the  Encyclopedia  cannot  provide 
all  the  necessary  details  on  the  problem.  Our  goal  is  to  provide  an  introduction  to  the  flow 
instability  concept  focusing  on  the  Tollmien-Schlichting  mechanism  in  boundary  layers  and 
to  outline  how  the  knowledge  about  this  mechanism  can  be  used  for  transition  prediction. 
Other  mechanisms  such  as  Gortler  and  cross-flow  mechanisms  are  discussed  in  in  the  books 
cited  above,  as  well  as  the  review  papers  [RS89,  Flo91,  Sar94a],  and  [SRW03]  on  these  topics. 
We  skip  the  discussion  of  free  shear  layers  (mixing  layers  and  wakes)  that  are  also  important 
in  technical  applications. 
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5.2  Stability  analysis  of  two-dimensional  flows 


5.2.1  Outline  of  the  linear  stability  theory 

For  simplicity,  consider  a  two-dimensional  flow  of  an  incompressible  fluid  governed  by  the 
Navier-Stokes  equations.  In  the  Cartesian  coordinates,  the  equations  are  written  as  follows: 


du  dv  ^ 
dx  dy 

du  du  du 

a + 

dv  dv  dv 

a 


1  dp  ( d'^u  d'^u  \ 

pdx  ^  \  dx"^  dy'^ ) 

1  dp  /  d'^v  d'^v  \ 

pdy~^  dy"^) 


(5.1) 


where  u  {x,  y,  t)  and  v  {x,  y,  t)  are  x—  and  y—  velocity  components,  respectively;  p  {x,  y,  t)  is 
the  pressure;  p  is  the  density;  and  z/  is  the  kinematic  viscosity.  We  assume  that  the  flow  field 
can  be  presented  as  the  sum  of  a  basic  steady  flow  and  a  perturbation: 


u  =  U  {x,y)+  u'  {x,  y,t),v  =  V  (x,  y)  +  v'  {x,  y,t),p  =  P  {x,  y)  +  p'  {x,  y,  t)  (5.2) 


After  substituting  the  decomposed  velocities  and  the  pressure  in  equations  (5.1)  and  sub¬ 
tracting  equations  for  the  basic  flow,  we  arrive  at  a  system  of  nonlinear  equations  for  per¬ 
turbations  m',  v' ,  and  p' .  In  the  linear  stability  analysis,  the  nonlinear  terms  are  neglected 
assuming  that  the  amplitude  of  the  perturbations  is  small. 
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(5.3) 


In  a  parallel  flow  such  as  a  fully  developed  channel  flow,  the  mean  velocity  profile  U  =  U  {y), 
V  =  0,  and  so  the  underlined  terms  in  equation  (5.3)  are  absent.  Boundary-layer  and 
free  shear  layers  (mixing  layer  and  far  wake)  are  not  parallel.  However,  these  flows  can  be 
considered  in  the  first  approximation  as  quasi-parallel  because  there  is  a  significant  difference 
in  the  length  scales.  There  is  a  length  scale,  L,  associated  with  the  weak  dependence  of  the 
mean  velocity  profiles  on  the  streamwise  coordinate  x.  Another  length  scale  is  the  boundary- 
layer  thickness  (thickness  of  the  mixing  layer  or  wake)  6.  At  high  Reynolds  numbers,  6/L 
I /Res  1,  where  Res  is  the  Reynolds  number  based  on  S  and  characteristic  velocity  (for 
example,  velocity  at  the  edge  of  the  boundary  layer,  Ue).  The  y—  velocity  component  of 
the  mean  flow,  V,  has  the  order  of  magnitude  Ue/Resi  and  U  ~  Ue-  In  addition,  the  length 
scale  of  the  perturbation.  A,  usually  has  the  same  order  of  magnitude  as  the  boundary-layer 
thickness  (L  3>  A  >  5).  As  a  result  of  the  scale  differences,  the  underlined  terms  in  equations 
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mode  1 


k 


mode  2  — ► 

Figure  5.1:  A  mechanical  system  having  two  degrees  of  freedom 


(5.3)  are  neglected  in  the  hrst  approximation  when  we  study  the  stability  of  boundary  and 
free  shear  layers.  In  this  case,  the  local  coordinate  x  in  the  mean  flow  velocity  prohle, 
U  {x,y),  is  treated  as  a  parameter.  In  other  words,  the  quasi-parallel  flow  approximation 
deals  with  the  local  velocity  prohles  that  are  considered  parallel  on  the  scale  having  an  order 
of  magnitude  A. 

The  quasi-parallel  flow  approximation  leaves  us  with  a  system  of  partial  (linear)  differen¬ 
tial  equations.  A  rigorous  mathematical  analysis  requires  the  formulation  of  the  initial  and 
boundary  conditions  for  the  perturbations  and  even  once  these  are  formulated,  the  prob¬ 
lem  may  remain  challenging.  In  practical  applications,  we  are  not  interested  in  the  solution 
of  specihc  initial- value  problems.  The  main  questions  of  interest  are:  1.  Is  the  flow  sta¬ 
ble/unstable?  2.  What  are  the  unstable  modes?  3.  Are  the  unstable  modes  amplihed  in  the 
flow  and  how  much? 

At  this  point,  we  should  address  the  normal  modes  concept.  The  simplest  illustration  of 
the  normal  modes  concept  is  the  mechanical  system  shown  in  hgure  5.1.  The  system  has  two 
masses  m  connected  by  a  spring  of  stiffness  k.  The  masses  are  also  attached  to  the  walls  by 
springs.  The  system  has  two  normal  modes  of  vibration:  Mode  1  has  frequency  u  = 
when  the  masses  oscillate  in  phase,  whereas  mode  2  has  frequency  on  =  (3/c/m) and  the 
displacements  of  the  masses  are  shifted  in  tt.  It  is  known  that  any  vibration  of  the  system 
can  be  presented  as  a  sum  of  these  normal  modes.  The  study  of  the  normal  modes  provides 
information  about  the  properties  of  the  system.  In  order  to  hnd  the  amplitudes  and  phases 
of  each  mode,  we  have  to  address  a  specihc  initial-value  problem  for  the  system. 

In  hydrodynamic  stability  analysis,  the  normal  mode  concept  suggests  considering  simple 
wave-like  solutions  in  the  complex  form  as  q'  {x,  y,t)  =  q  {y)  exp  {iax  —  iut),  where  q  {y)  is 
a  complex  amplitude  function.  The  idea  is  the  same  as  in  the  example  of  the  mechanical 
system:  instead  of  solving  a  specihc  physical  initial  and  boundary  value  problem,  we  look 
at  the  normal  modes  existing  in  the  system,  expecting  that  the  solution  of  the  initial  and 
boundary  value  problem  for  the  partial  diherential  equations  (5.3)  can  be  presented  as  a  sum 
of  the  normal  modes.  If  there  is  an  unstable  mode,  one  can  expect  that  the  mode  will  also 
be  present  in  the  solution  of  a  specihc  physical  problem  (in  a  real  experiment  or  in  direct 
numerical  simulation)  and  that  the  mode  can  be  observed  (perhaps  after  a  space  or  time  delay 
required  for  sufficient  amplihcation  of  the  mode).  The  main  diherence  between  the  simple 
mechanical  system  and  perturbations  in  a  huid  is  that  the  continuum  medium  has  an  inhnite 
number  of  degrees  of  freedom.  As  result,  one  can  arrive  at  a  enumerable  discrete  spectrum 
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or  at  a  continuous  spectrum,  or  at  a  combination  of  discrete  and  continuous  spectra. 

Following  the  normal  modes  concept,  the  solution  of  equations  (5.3)  is  sought  in  the  form 

(m',  n',  p')  =  (h,  h,  p)  exp  {iax  —  iut)  (5.4) 


and  the  system  of  the  partial  differential  equations  is  reduced  to  a  system  of  ordinary  differ¬ 
ential  equations  (ODEs)  for  the  amplitude  functions. 
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Although  we  use  the  symbol  of  an  ordinary  derivative  dU / dy  in  the  second  equation,  in  the 
case  of  weakly  non-parallel  flows  we  use  the  velocity  prohle  at  a  chosen  coordinate  x  that  is 
considered  as  a  parameter. 

The  normal  modes  are  described  by  the  complex  amplitude  functions  that  provide  am¬ 
plitudes  and  phases  depending  on  the  coordinate  y.  Similarly  to  the  spring-mass  problem 
discussed  above,  solution  of  the  linearized  Navier-Stokes  equations  can  be  presented  as  an 
expansion  into  the  normal  modes. 

The  system  of  equations  (5.5)  can  be  recast  as  one  ODE  of  the  fourth  order  for  the 
amplitude  function  v: 


{all  —  u)  {y”  —  a^v)  —  avU"  =  —iv  (h™  —  2a^v"  -|-  a^v)  (5.6) 

where  (...)^  stands  for  d/dy,  and  =  d^v/dy^.  Usually,  a  perturbation  of  streamfunction 
(p  iy)  exp  {iax  —  iuit)  is  introduced  {v  =  —iap),  and  the  equation  (5.6)  is  written  in  dimen¬ 
sionless  form  (using  the  characteristic  length  scale  5  and  the  velocity  scale  Ue)  as  follows: 

{alJ  —  u)  [if"  —  a^(f)  —  acpU"  =  — —  —  2a‘^(p"  +  a^^cp^  (5.7) 

This  is  the  famous  Orr-Sommerfeld  equation  in  hydrodynamic  stability  theory.  One  may 
compare  the  Orr-Sommerfeld  equation  with  the  system  of  ODEs  (5.5)  written  using  the 
primitive  variables  (velocity  components  and  pressure).  Equation  (5.7)  requires  the  second 
derivative  of  the  mean  flow  prohle,  whereas  equations  (5.5)  contain  only  the  hrst  deriva¬ 
tive.  Because  theoretical  studies  historically  had  been  linked  to  equation  (5.7),  there  is  a 
misunderstanding  that  stability  analysis  of  low  speed  boundary  layers  requires  the  second 
derivative  U". 

Equation  (5.7)  (or  equations  (5.5)  )  require  boundary  conditions.  On  a  solid  wall,  the 
no-slip  boundary  condition  must  be  satished: 


u  =  p'  =  0,  V  =  —iap  =  0 


(5.8) 
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In  the  case  of  an  nnbounded  domain  (as  in  boundary-layer  flow),  one  may  require  bound¬ 
edness  of  the  solution  outside  the  boundary  layer.  One  can  hnd  a  solution  that  decays 
outside  the  boundary  layer  {u,v,p  — >  0),  and  there  is  also  the  option  that  the  solution  is 
bounded  but  does  not  decay  (having  an  oscillating  character  along  coordinate  y).  The  hrst 
case  corresponds  to  the  discrete  spectrum,  and  the  second  option  leads  to  the  continuous 
spectrum.  The  solutions  from  the  continuous  spectrum  can  be  associated  with  perturbations 
in  the  free  stream.  For  example,  turbulence  in  the  free  stream  may  be  associated  with  a 
normal  mode  of  the  continuous  spectrum.  In  the  case  of  a  compressible  gas,  acoustic  and 
entropy  perturbations  in  the  free  stream  can  be  presented  by  the  continuous  spectrum. 

In  what  follows,  our  main  interest  will  be  associated  with  modes  of  the  discrete  spectrum 
because  these  include  the  unstable  modes.  In  order  to  clarify  the  origin  of  the  discrete 
spectrum,  consider  equation  (5.7)  outside  the  boundary  layer,  where  the  dimensionless  mean 
flow  velocity  U  =  1.  At  y  ^  oo,  the  Orr-Sommerfeld  equation  becomes  an  equation  with 
constant  coefficients 

{a  —  u)  {ip"  —  a^ip')  =  — —  {ip™  —  ‘lo?ip"  -|-  a^ip')  (5.9) 

6'XLC'(5 

and  one  can  easily  hnd  four  linear  independent  solutions 

ipi  =  exp  {-ay) ,  ip2  =  exp  {+ay) ,  993  =  exp  (-7?/) ,  ip^  =  exp  (-^7?/)  (5.10) 

7  =  -1-  iRes  {a  —  u) 

To  be  specihc,  we  chose  real  parts  of  a  and  A  to  be  positive.  Therefore,  only  two  linear 
independent  solutions  ipi  and  993  decay  outside  the  boundary  layer.  The  decaying  solution 
can  be  written  as  a  sum 

ip{y)  =  Cip>i{y)  +  C^p>^{y)  (5.11) 

where  constants  Ci  and  are  unknown.  One  can  continue  ipi  and  993  into  the  boundary 
layer  numerically  and  try  to  hnd  the  constants  Ci  and  Os  from  the  wall  boundary  conditions 
(5.8).  However,  one  of  the  constants  can  be  chosen  arbitrarily  because  the  equation  is  linear. 
The  choice  of  the  constant  will  ahect  the  normalization  of  solution  (5.11)  only.  The  other 
constant  can  be  found  from  a  boundary  condition  on  the  wall,  whereas  the  second  wall¬ 
boundary  condition  still  cannot  be  satished.  Only  a  special  choice  of  parameters  Q;,a;,  and 
Res  can  lead  to  the  solution  satisfying  both  boundary  conditions  (5.8).  This  is  the  discrete 
spectrum  obeying  a  dispersion  relation: 

uj  =  u{a,Res)  OTa  =  a{uj,Res)  (5-12) 

At  this  point,  the  normal  mode  analysis  leads  to  an  ambiguity  that  can  be  readily  found  in 
the  literature.  In  the  general  case,  we  may  consider  the  wave  number,  a  =  ar  +  iai,  and 
frequency,  u  =  Ur  +  iooi  as  complex  numbers.  Very  often  two  formulations  are  mentioned: 
temporal  and  spatial  theories. 

In  the  temporal  stability  theory,  Oj  =  0  and  is  considered  as  a  parameter,  while  the 
complex  number  u  is  found  from  the  dispersion  relation  (5.12).  According  to  the  normal 
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mode  presentation,  the  solution  depends  on  time  as  exp(— icut)  =  exp(— +  Uit).  At 
Ui  >  0,  the  solution  is  unstable  (it  grows  exponentially  in  time),  and  it  decays  at  ct;*  <  0  (the 
solution  is  neutral  at  Ui  =  0). 

In  the  spatial  theory,  Ui  =  0  and  Ur  is  considered  as  a  parameter,  while  the  complex 
number  a  is  found  from  the  dispersion  relation  (5.12).  Now,  the  amplitude  of  the  solution 
depends  on  the  coordinate  x  as  exp(iQ;a;)  =  exp{iarX  —  aix) .  At  a*  <  0,  the  amplitude  grows 
exponentially  downstream,  and  the  perturbation  decays  downstream  at  >  0. 

The  reader  can  hnd  an  ambiguity  in  the  literature  within  the  choice  of  temporal  or  spatial 
formulation  in  the  stability  problems.  Usually,  the  spatial  formulation  is  recommended  as 
more  suitable  for  analysis  of  experimental  data.  The  ambiguity  has  its  origin  at  the  step 
when  the  solution  of  the  partial  differential  equations  is  suggested  in  the  form  of  the  normal 
modes  (5.4).  At  this  point,  the  link  with  the  initial  and  boundary  value  problem  for  a  specihc 
experimental  setup  is  lost,  and  we  have  the  artihcial  dilemma  of  spatial  versus  temporal 
formulations. 

[Gus79]  solved  an  initial  value  problem  for  a  perturbation  in  an  incompressible  boundary 
layer.  [SG81]  proved  that  Gustavsson’s  solution  is  an  expansion  into  the  modes  of  discrete 
and  continuous  spectra  in  the  temporal  formulation.  [AR90]  considered  the  signaling  problem 
when  a  periodic-in-time  point  source  introduces  perturbations  through  the  wall  starting  at 
t  =  0.  They  showed  that  the  solution  of  the  initial  and  boundary  value  problem  can  be 
presented  (after  some  time  when  the  transient  effect  is  smeared  out)  as  an  expansion  into 
the  modes  of  discrete  and  continuous  spectra  in  the  spatial  formulation.  Therefore,  the 
analysis  of  a  specihc  initial  boundary  value  problem  for  PDEs  leads  to  a  unique  expansion 
of  the  solution  into  normal  modes  without  any  ambiguity.  Usually  experiments  deal  with 
disturbance  source  like  in  the  signaling  problem,  and  the  spatial  framework  should  be  used 
for  analysis  of  the  data. 

5.2.2  Blasius  boundary- layer  flow 

The  neutral  stability  curve  for  a  Blasius  boundary  layer  how  (section  refeae012)  is  shown  in 
hgure  5.2,  where  F  =  (/  is  the  dimensional  frequency)  and  Re  =  {Uex/i')^^'^ . 

Figure  5.2  illustrates  the  remarkable  two-fold  role  of  viscosity.  The  how  is  stable  at  low 
Reynolds  numbers  due  to  the  viscous  dissipation.  The  how  is  also  stable  at  large  Reynolds 
numbers.  This  means  that  the  underlying  instability  mechanism  is  viscous  (see  a  discussion 
in  section  5.2.4).  The  unstable  modes  having  a  viscous  nature  are  called  Tollmien-Schlichting 
(TS)  waves  to  honor  the  theoretical  studies  of  W.  Tollmien  and  H.  Schlichting  in  the  late 
1920s  and  the  early  1930s  (see  the  references  in  [SGOO]).  Their  theoretical  predictions  were 
conhrmed  by  the  experiments  of  [SS43].  Stability  experiments  require  quiet  wind  tunnels  and 
high  accuracy  of  measurements.  The  reader  can  hnd  a  number  of  useful  data  and  description 
of  recent  experiments  in  review  papers  [Sar94b]  and  [Sar07]. 

The  imaginary  part  of  the  wave  numbers  (in  the  spatial  formulation)  are  shown  in  hg¬ 
ure  5.3  as  functions  of  the  Reynolds  number  at  the  frequency  parameters  F  =  n  x  10“^ 
(n  =  1,. .  .,5).  The  Blasius  length  scale  Lb  =  (ra/Ug)^'^^  is  used  in  the  dehnition  of  the 
dimensionless  wave  number. 
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Figure  5.2:  Neutral  stability  curve  for  Blasius  boundary  layer  flow  in  the  quasi-parallel  flow 
approximation. 


Figure  5.3:  Imaginary  part  of  the  wave  numbers  versus  the  Reynolds  number  for  Blasius 
boundary  layer  flow.  The  numbers  n  =  1,  2,  3,  4,  and  5  indicate  the  frequency  parameters 
F  =  n  X  10“^. 
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Figure  5.4:  Amplitude  distributions  for  streamwise  velocity  perturbation  at  F  =  5  x  10“® 


A  typical  streamwise  velocity  profile  of  the  TS  waves  in  a  Blasius  boundary  layer  is  shown 
in  hgure  5.4  at  F  =  5  x  10“®  at  two  Reynolds  numbers,  where  the  dimensionless  distance 
from  the  wall  is  dehned  using  the  Blasius  length  scale,  rj  =  y/ Lb-  A  comparison  of  the 
velocity  prohles  in  hgure  5.4  illustrates  the  non-parallel  how  ehect:  the  shape  of  the  velocity 
prohles  depends  on  the  distance  from  the  leading  edge.  One  can  consider  an  experiment 
where  velocity  perturbations  are  measured  using  a  hot  wire.  One  can  see  from  hgure  5.4 
that  the  measured  growth  (decay)  of  the  perturbation  will  depend  on  the  location  of  the 
hot  wire  within  the  boundary  layer.  In  addition  to  quantitative  corrections  required  when 
the  nonparallel  how  ehects  are  taken  into  account,  the  perturbation  growth  rate  will  depend 
on  the  measured  quantity.  For  example,  one  may  consider  amplihcation  measured  for  the 
inner  or  outer  maxima  in  the  velocity  prohle.  The  nonparallel  how  ehects  in  incompressible 
hows  were  studied  using  the  method  of  multiple  scales  by  [Bou72,  Gas74];  [SN75,  SN77];  and 
[PN79].  Figure  5.5  illustrates  the  neutral  curves  ([Gas74])  dehned  for  the  inner  and  outer 
maxima  of  the  streamwise  velocity  prohle  in  a  TS  wave.  The  Reynolds  number  in  hgure 
5.5  is  dehned  using  the  displacement  thickness  5* .  Experimental  studies  of  non-parallel  how 
ehects  are  discussed  by  [Sar07]. 

5.2.3  The  Falkner-Skan  profiles.  The  pressure  gradient  effect 

The  Falkner-Skan  boundary  layer  how  is  characterized  by  the  edge  velocity  Ue{x)  =  Cx^. 
In  this  case,  it  is  common  to  use  the  Hartree  parameter  (3  =  2m/{m  +  1)  (section  refeae012). 
Stability  analyses  have  shown  that  a  favorable  pressure  gradient  {(3  >  0)  stabilizes,  and  an 
unfavorable  pressure  gradient  (/?  <  0)  destabilizes  the  boundary  layer  how.  [DR81]  collected 
data  for  the  critical  Reynolds  numbers  (minimal  Reynolds  numbers  on  the  neutral  curve) 
of  Falkner-Skan  prohles  from  [OML69]  and  plotted  them  as  function  of  the  shape  factor 
H  =  5* / 9,  where  9  is  the  momentum  thickness.  Figure  5.6  reproduces  the  plot  from  [DR81]. 
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Figure  5.5:  Neutral  stability  curves  for  Blasius  boundary  layer  flow  (Gaster,  1974.  Repro¬ 
duced  by  permission  of  Cambridge  University  Press).  Solid  line  quasi-parallel  flow  approx¬ 
imation;  dash  line  inner  maximum  of  |u|;  dash-dot  line  outer  maximum  of  |m|. 


Figure  5.6:  The  values  of  the  minimal  critical  Reynolds  number  (based  on  the  displacement 
thickness)  for  the  Falkner-Skan  family  of  boundary  profiles  (Drazin  and  Reid,  1981). 
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5.2.4  The  Rayleigh  equation  and  structure  of  the  TS  wave  at  high 
Reynolds  numbers 

In  the  limit  of  high  Reynolds  numbers,  the  solution  of  the  Orr-Sommerfeld  equation  can  be 
considered  in  a  rigorous  way  using  the  method  of  matched  asymptotic  expansions  ([Smi79, 
SB80,  BS81,  ZR83]).  In  the  following  discussion,  we  are  using  a  simplihed  structure  of  the 
solution  at  high  Reynolds  numbers  in  order  to  outline  the  main  ideas  ([Gra66,  Eag69,  DR81, 
SHOl]). 

In  the  limit  aRes  oo,  one  can  ignore  the  viscous  terms  in  equation  (5.7).  The  inviscid 
approximation  is  the  Rayleigh  equation: 

{U  -  c)  {(p"  -  a^(p)  -  (pU"  =  0  (5.13) 

where  c  =  u/ a  is  the  phase  velocity.  One  can  expect  that  the  viscous  terms  are  still  signihcant 
in  the  vicinity  of  the  wall  (viscous  sublayer)  in  order  to  satisfy  the  no-slip  boundary  condition. 
In  addition,  equation  (5.13)  is  singular  in  the  vicinity  of  the  point  i/c,  where  U  (i/c)  =  c  (critical 
layer).  The  viscous  terms  are  significant  in  the  vicinity  of  the  critical  layer  as  well.  In  order 
to  hnd  an  estimate  of  the  viscous  sublayer  thickness  S^s,  we  compare  the  order  of  magnitude 
of  the  viscous  term  having  the  highest  derivative  with  the  order  of  magnitude  of  the  term 
{U  —  c)(p"  at  y  — 0: 


{U  -  c)ip"  ~  cp>”  ~  c(p/5l^, 


aRes 


(p  ~ 


aResS^g 


(5.14) 


These  terms  are  balanced  when  ~  (acRes)  In  order  to  estimate  the  thickness  of  the 
critical  layer,  we  consider  a  balance  of  the  same  terms  at  y  ^  yc- 


{U  -  c)(p”  ~ 


U',p>5, 

SI 


aRes 


aRes5* 


(5.15) 


where  U'c  is  the  derivative  at  the  critical  point.  A  balance  of  the  terms  requires  5c  ~ 
[alJ' cRgs)~^^^ ■  There  are  two  possibilities.  The  critical  layer  can  be  close  to  the  wall  and 
overlap  with  the  viscous  sublayer.  In  this  case,  the  viscous  effects  are  signihcant  in  the 
vicinity  of  the  wall  only.  Such  a  structure  corresponds  to  the  asymptotic  behavior  of  the 
lower  branch  of  the  neutral  curve.  Another  possibility  is  when  the  critical  layer  and  viscous 
sublayer  are  separated,  and  there  is  an  inviscid  region  between  the  viscous  layers.  This 
structure  corresponds  to  the  upper  branch  of  the  neutral  curve. 

In  order  to  shed  some  light  onto  the  nature  of  the  viscous  instability  mechanism,  let  us 
consider  a  neutral  perturbation  with  a  >  0  corresponding  to  the  upper  branch  of  the  neutral 
curve.  One  can  easily  derive  from  the  Rayleigh  equation  the  following  result: 

^  {cpcp' -  cpcp')  =  0  (5.16) 

{(p(p'  —  pxp')  =  2i  Im  {<p(p')  =  const 
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where  bar  stands  for  complex  conjugation. 

In  order  to  understand  the  physical  meaning  of  {(pip'  —  ipip'),  consider  the  Reynolds  stress 
that  the  neutral  perturbation  generates  in  the  boundary  layer. 

T 

T  =  -^j  u*  {x,y,t)v*  {x,y,t)dt  (5.17) 

0 

T  =  27r/ui 

The  velocity  components  u*  {x,  y,  t)  and  v*  (x,  y,  t)  are  presented  using  the  complex  stream 
function  as  follows 


u*  {x,  y,  t)  =  Re  ,  v*  {x,  y,  t)  =  Re  (5.18) 

After  substitution  of  (5.18)  into  integral  (5.17),  one  can  find  the  Reynolds  stress 

^  [w' -  Wl  =  |lm((^(^')  (5.19) 


Using  equation  (5.16),  we  conclude  that  the  Reynolds  stress  is  constant  within  the  inviscid 
regions  at  aRes  oo.  Because  the  eigenfunction  corresponding  to  the  TS  wave  decays  at 
y  — cx),  we  conclude  that  the  Reynolds  stress  r  =  0  in  the  inviscid  region  above  the  critical 
layer.  The  Reynolds  stress  also  has  to  be  constant  within  the  inviscid  region  between  the 
viscous  sublayer  and  the  critical  layer,  and  the  sign  of  r  is  very  important.  In  order  to  have 
energy  production,  r  must  be  positive. 

We  can  explore  the  sign  of  r  using  an  approximate  solution  of  the  Orr-Sommerfeld  equa¬ 
tion  in  the  vicinity  of  the  wall.  The  leading  term  is  determined  by  the  solution  of  the 
following  equation: 


-uip"  = 


iRes 


(5.20) 


Equation  (5.20)  was  obtained  from  equation  (5.7)  when  the  highest  derivatives  are  kept  in 
viscous  and  inviscid  terms  aX  y  ^  0.  The  no-slip  boundary  conditions  on  the  wall  are  as 
follows: 


43(0)=43'(0)  =  0 


(5.21) 


The  solution  of  the  equation  satisfying  the  boundary  conditions  can  be  easily  found  (linear 
independent  solution  exp{+q()  has  to  be  excluded): 


ip  =  e  +  q(  —  1 
1  —  i 


,  C  =  y{uResY^‘^ 


In  the  limit  (  ^  ryo  (outer  limit  of  the  inner  solution)  we  hnd. 


(5.22) 


ip^  qC-1 


(5.23) 
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Figure  5.7:  Sketch  of  the  Reynolds  stress  distribution  across  the  boundary  layer  in  the  high 
Reynolds  number  limit  when  the  viscous  sublayer  and  critical  layer  are  separated  (upper 
branch  of  the  neutral  curve).  1  viscous  sublayer;  2  inviscid  region;  3  critical  layer;  4 
inviscid  region. 


Substitution  of  the  result  into  (5.19)  gives  the  Reynolds  stress  outside  of  the  viscous  sublayer. 


ayJujRe^ 

2^2 


(5.24) 


A  sketch  of  the  asymptotic  Reynolds  stress  structure  is  shown  in  hgure  5.7.  It  is  important 
that  r  is  positive.  In  other  words,  the  viscous  sublayer  provides  a  phase  shift  between 
the  streamwise,  u*  {x,y,t),  and  normal,  v*  {x,y,t),  velocity  components  leading  to  r  >  0 
(energy  production).  In  the  case  of  the  neutral  perturbation,  one  can  show  that  the  energy 
production  between  the  critical  and  viscous  sublayer  is  balanced  by  the  energy  dissipation 
within  the  viscous  sublayer  ([ZT87]). 


5.2.5  The  Rayleigh  theorem 

Using  the  Rayleigh  equation,  one  can  establish  the  following  Rayleigh  theorem  ([DR81, 
SHOl]):  A  necessary  condition  for  inviscidly  unstable  flow  is  the  existence  of  an  inflection 
point  U"  {y^,)  =0.  The  Rayleigh  theorem  provides  a  necessary  but  not  a  sufficient  condition. 
It  is  interesting  to  compare  neutral  curves  for  Blasius  boundary  layer  flow  with  neutral  curves 
for  boundary  layers  with  an  unfavorable  pressure  gradient  when  the  velocity  prohle  has  an 
inflection  point.  A  sketch  of  the  neutral  curves  is  shown  in  hgure  5.8  ([Pan05]).  Other  typical 
examples  of  inviscidly  unstable  hows  are  mixing  layers  and  wakes.  Examples  of  their  stability 
analysis  are  presented  in  books  by  [SHOl],  and  [CJJ67].  Although  U"  is  not  required  for 
the  stability  analysis  (section  5.2),  an  inhection  point  indicates  a  possibility  for  the  inviscid 
mechanism  of  the  how  instability. 

Inviscid  stability  analysis  of  a  mixing  layer  becomes  very  simple  in  the  case  when  the 
wave  length  of  the  perturbation  is  much  larger  than  the  mixing  layer  thickness  (A  3>  5). 
In  this  limit,  one  may  consider  two  uniform  streams  having  velocities  Ui  and  U2-  The 
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Figure  5.8:  General  shape  of  the  neutral  stability  curve  for  a  Blasius  boundary  layer  and  for 
layers  with  adverse  pressure  gradient  (Panton,  2005). 

slip  surface  is  the  sheet  vortex  of  uniform  density.  The  stability  analysis  includes  hnding 
perturbations  of  the  velocity  potential  in  the  uniform  streams.  The  pressure  perturbations 
and  normal  velocities  of  particles  at  the  disturbed  interface  are  to  be  equal  on  both  sides. 
These  conditions  lead  to  a  simple  result  for  the  complex  phase  velocity 

c=\{Ui  +  U2)±'-\U2-U^  (5.25) 

Sign  “+”  in  equation  (5.25)  corresponds  to  unstable  mode.  The  result  means  that  sheet 
vortex  is  unstable.  This  is  the  Kelvin-Helmholtz  instability  (see  [Bat03,  Pan05]).  Compu¬ 
tational  studies  of  mixing  layers  demonstrate  that  the  hnite  thickness  5  ~  A  stabilizes  the 
flow.  However,  the  flow  is  always  unstable  in  the  limit  A  3>  5  when  Ui  ^  U2- 

5.2.6  Three-dimensional  perturbations 

The  normal  mode  analysis  in  the  case  of  three-dimensional  perturbations  includes  depen¬ 
dence  on  the  spanwise  coordinate  ;2.  Three  velocity  components  and  pressure  perturbations 
are  sought  in  the  form 

(«',  v',  w',  p')  =  {it,  V,  w,  p)  exp  {iax  +  i^z  —  iuit)  (5.26) 

After  substitution  of  (5.26)  into  linearized  Navier-Stokes  equations,  one  can  derive  a  system 
of  ODEs  for  the  amplitude  functions  {u,v,w,p)  similar  to  equation  (5.5).  One  can  also 
derive  one  ODE  of  the  fourth  order  for  the  amplitude  function  v: 

{U  —  c)  (v”  —  -I-  /5^)  h)  —  vU"  =  - fh™  —  2  -|-  v"  +  (5.27) 

[Squ33]  found  that  one  can  transform  the  equation  into  the  form  of  an  Orr-Sommerfeld 
equation  for  two-dimensional  disturbances: 

Oi2D  =  Oi2D^^&2D  =  (5.28) 
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Transformation  (5.28)  is  the  Sqnire  transformation  leading  to  the  two-dimensional  form  of 
the  Orr-Sommerfeld  eqnation: 


(U 


c  [V 


0120^ 


))  -  W"  = 


ia2D^^&2D 


v"  + 


2D 


(5.29) 


Using  the  temporal  framework  of  the  stability  analysis,  one  can  say  the  following  abont 
three-dimensional  pertnrbations  ([Squ33]):  “if  any  velocity  prohle  is  unstable  for  a  partic¬ 
ular  value  of  Reynolds’  number,  it  will  be  unstable  at  a  lower  value  of  Reynolds’  number 
for  two-dimensional  disturbances.”  Squire’s  transformation  does  not  have  a  straight-forward 
spatial  formulation;  it  leads  to  complex  values  of  Re52D  that  defy  simple  physical  inter¬ 
pretation.  And,  in  fact,  one  can  hnd  (at  low  frequencies)  that  the  lower  branch  of  the 
neutral  curve  for  three-dimensional  perturbations  can  be  be  shifted  upstream  with  respect 
to  the  neutral  curve  for  two-dimensional  perturbations.  However,  the  overall  amplihcation 
of  two-dimensional  disturbances  is  larger  than  the  amplihcation  of  the  three-dimensional 
perturbations  having  the  same  frequency.  Although  the  fact  is  often  interpreted  as  a  re¬ 
sult  of  the  Squire  transformation,  this  is  not  quite  correct.  In  contrast,  three-dimensional 
disturbances  are  more  amplihed  in  supersonic  boundary  layers. 

5.3  Receptivity  of  boundary  layers 

Having  introduced  the  normal  modes  concept  in  section  5.2.1,  we  pointed  out  that  the 
analysis  of  the  normal  modes  is  missing  a  link  with  the  PDEs  and  the  physical  initial  and 
boundary  value  problem.  In  other  words,  the  normal  mode  analysis  substitutes  solving  a 
real  physical  problem  with  studies  of  some  elements  of  the  bases  (eigenfunctions),  omitting  a 
discussion  of  what  their  weights  might  be  in  a  physical  setup.  (This  problem  is  often  called 
the  receptivity  problem).  It  had  been  understood  by  the  mid-1970s  that  the  instability  of 
flow  is  only  one  element  of  the  laminar-turbulent  transition  scenario.  [Mor69]  and  [Res76] 
clarihed  the  important  role  of  the  receptivity  problem  in  the  laminar-turbulent  transition 
process.  These  papers  motivated  intensive  investigations  of  various  mechanisms  responsible 
for  the  excitation  of  unstable  TS  waves.  The  vast  bibliography  on  the  topic  is  presented 
in  [Cho98,  SRK02,  CJJ67,  Fed03];  and  [TumOOb,  Tum07].  For  our  further  discussion  it  is 
important  to  keep  in  mind  that  the  receptivity  issue  is  a  signihcant  element  in  the  formulation 
of  a  physics-based  transition  prediction. 


5.4  Paths  to  turbulence  in  wall  layers 

Before  a  discussion  of  possible  scenarios  for  laminar-turbulent  transition,  we  have  to  point 
out  that  the  previous  sections  were  aimed  at  modes  that  can  be  unstable  and  lead  to  the 
transition.  However,  this  concept  is  not  adequate  for  the  whole  spectrum  of  experimental 
observations.  For  example,  the  Couette  flow  and  Hagen-Poiseuille  flow  in  a  circular  pipe  are 
stable,  though  one  can  observe  in  experiments  a  laminar-turbulent  transition  in  these  flows. 
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There  are  also  examples  when  a  transition  to  turbulent  flow  can  be  triggered  at  Reynolds 
numbers  lower  than  the  critical  value  predicted  by  the  stability  theory.  These  observations 
indicate  that  there  are  other  routes  in  transition  to  turbulence  aside  from  amplification  of  the 
instability  modes.  One  of  the  possible  mechanisms  of  transition  to  turbulence  is  associated 
with  the  transient  growth  phenomenon. 

Small  perturbations  in  a  stable  shear  flow  can  have  a  significant  transient  amplification 
before  they  die  out  in  accordance  with  the  stability  theory  prediction.  A  pair  of  counter¬ 
rotating  streamwise  vortices  in  a  boundary  layer  illustrates  the  transient  growth  phenomenon. 
The  vortices  redistribute  the  streamwise  momentum  by  lifting  up  slow  fluid  particles  from  the 
wall  and  bringing  down  the  high-speed  particles  toward  the  wall.  Although  the  streamwise 
vortices  decay  downstream,  the  cumulative  effect  of  the  momentum  redistribution  can  be 
significant,  and  one  can  observe  the  streaky  structures  in  boundary  layers.  In  turn,  the 
transiently  growing  perturbations  can  lead  to  a  flow  structure  that  is  unstable  with  respect 
to  other  small  perturbations,  serving  as  a  path  to  turbulence  without  the  presence  of  unstable 
modes  (non-modal  path  to  transition).  One  can  And  bibliographies  on  theoretical  studies  of 
this  phenomenon  in  [SHOl],  and  [ZTR06].  The  transient  growth  phenomenon  can  be  observed 
in  experiments  with  an  array  of  roughness  elements  placed  on  the  wall  ([Whi02,  FBTC04, 
WRE05,  EW06,  DW08]).  The  transient  growth  phenomenon  has  been  studied  extensively 
using  theoretical,  experimental,  and  computational  methods  in  recent  years  because  of  its 
potential  application  for  transition  prediction  ([ResOl,  RT04,  RT06]. 

Professor  E.  Reshotko  ([RT06])  analyzed  possible  paths  in  transition  to  turbulence  and 
summarized  the  contemporary  understanding  in  the  diagram  shown  in  figure  5.9.  It  is  impor¬ 
tant  to  point  out  that  the  path  to  turbulence  depends  on  the  level  of  external  disturbances. 
Depending  on  the  environment  (noisy  wind  tunnel,  flight  in  a  quiet  atmosphere,  flow  in  a 
turbine  etc),  different  paths  can  have  different  values  at  particular  conditions,  and  a  physics- 
based  transition  prediction  has  to  take  into  account  the  variety  of  mechanisms. 

Despite  the  variety  of  the  transition  mechanisms,  they  all  lead  to  nonlinear  flow  dynamics 
at  the  final  stage  of  the  flow  breakdown.  Historically,  the  following  transition  scenarios  have 
been  distinguished:  the  K- regime  (after  [KTS62]),  the  subharmonic-  (also  known  as  N-  or 
H-)  regime  discovered  by  [KKL77],  and  the  oblique  transition  studied  by  [FTB93,  BLH94], 
and  others.  Discussions  of  the  nonlinear  mechanisms  in  transition  scenarios  can  be  found 
in  [Kac94],  and  [SHOl].  Recent  experimental  and  computational  studies  indicate  that  these 
different  mechanisms  are  the  same  at  the  late  stages  of  the  transition  ([SFKOO]).  Experiments 
in  a  circular  pipe  flow  ([HTWOO])  also  demonstrate  that  the  late  stage  of  transition  is  similar 
to  transition  in  boundary  layers. 


5.5 


Transition  prediction:  modal  growth  scenario  (path 

A) 


The  laminar-turbulent  transition  depends  on  many  factors,  and  prediction  of  the  “transition 
point”  is  a  challenging  task.  Because  the  phenomenon  depends  on  external  forcing  (free 
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increasing  disturbance  level 


Figure  5.9: 


Paths  to  turbulence  in  wall  layers  (Reshotko  and  Tumin,  2006). 


stream  turbulence,  acoustic  noise,  roughness  induced  perturbations  etc),  transition  predic¬ 
tion  suffers  uncertainty  related  to  the  ffow  conditions.  The  level  of  perturbations  (and  their 
spectrum)  in  wind  tunnels  might  be  different  from  flight  conditions,  and  one  cannot  use  wind 
tunnel  experiments  for  transition  prediction  in  a  flight.  However,  if  the  transition  path  is  de¬ 
termined,  one  can  work  on  a  physics-based  transition  prediction  method.  Consider  transition 
on  a  smooth  snrface  dne  to  low-level  free  stream  pertnrbations  (path  A  in  fignre  5.9).  The 
main  stages  of  the  transition  are:  receptivity,  linear  amplification  of  TS  waves,  and  nonlinear 
breakdown.  As  follows  from  experiments  in  qniet  wind  tunnels  and  their  comparisons  with 
the  linear  stability  theory,  a  significant  part  of  the  transition  region  can  be  described  by  the 
linear  stability  theory,  and  the  length  of  the  nonlinear  region  is  relatively  short.  Comparison 
of  experimental  data  with  the  theoretical  results  (see,  for  example,  [ZT87])  indicates  that  the 
nonlinear  dynamics  of  perturbations  become  significant  when  the  velocity  of  the  perturba¬ 
tion  achieves  abont  1%  of  the  free  stream  velocity.  Therefore,  one  can  assnme  that  there  is  a 
threshold  amplitude  e*  associated  with  the  laminar-turbulent  transition  (amplitude  criterion 
of  transition).  Let’s  assnme  that  we  know  the  spectrnm  of  the  free  stream  pertnrbations, 
£oo(/),  and  the  receptivity  coefficient  K{f)  snch  that  the  amplitnde  of  the  nnstable  wave  at 
the  nentral  point,  Xnp,  can  be  fonnd  as  ersif)  =  -^(/)^oo(/)-  The  linear  amplification  of 
the  pertnrbation  can  be  predicted  by  the  factor  exp(—  Ui  dx).  Therefore,  the  “transition 
point,”  Xtr  can  be  estimated  from  the  amplitude  criterion: 

=  K{f)eooif)exp{-  aidx)  (5.30) 

^np 

Equation  (5.30)  can  be  recast  as  follows: 
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Figure  5.10:  iV-factors  for  the  Blasius  boundary  layer  flow  (quasi-parallel  flow  approxima¬ 
tion). 


N{xtr) 

N{xtr) 


In 


K{f)eoo{f) 

r^tr 

I  ai  dx 

^np 


(5.31) 


Because  the  effects  of  the  unknown  factors,  K{f)eoo{f),  are  smoothed  by  the  logarithm,  one 
can  try  to  correlate  N{xtr)  for  different  flow  conditions  (noisy/quiet  wind  tunnels,  flight, 
etc).  This  is  the  underlying  concept  of  the  method  transition  prediction  ([vI56,  Smi59]). 
Figure  5.10  illustrates  the  iV-factors  for  the  Blasius  boundary  layer  in  the  quasi-parallel 
approximation. 

Because  the  nonparallel  flow  effect  can  have  a  signihcant  impact  on  the  growth  rate, 
the  contemporary  applications  of  the  method  are  based  on  the  Parabolized  Stability 
Equations  (PSE)  introduced  by  [HB87]  (see  also  [BerOl,  Her97]).  In  the  PSE  method,  one 
has  to  solve  a  system  of  parabolized  equations  using  a  downstream  marching  algorithm.  The 
approach  allows  capturing  the  nonparallel  flow  effects  in  the  leading  term  without  additional 
steps  as  required  in  the  method  of  multiple  scales.  The  reader  can  End  a  bibliography  and 
examples  on  the  method’s  application  in  [HLRA96,  JC99,  JC05];[AC00,  SRW03,  CCNOl]; 
[Cha03,  Cha04];  and  [CC05]. 


5.6  Conclusion 

This  brief  discussion  demonstrates  that  transition  to  turbulence  is  a  complex  phenomenon. 
It  depends  on  a  variety  of  factors  (free-stream  disturbances,  perturbations  originated  on  the 
wall  etc)  that  might  be  beyond  our  control.  Therefore,  the  designer  of  laminarized  aircrafts 
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has  to  deal  with  the  uncertainties  and  to  distinguish  different  mechanisms  responsible  for 
the  flow  transition. 
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Chapter  6 


Toward  the  foundation  of  a  global 
(bi-global)  modes  concept 

This  section  is  based  on  the  invited  talk  at  the  fourth  symposium  on  Global  Flow  Instability 
and  Control  [Tum09]. 

6.1  Introduction 

Recent  progress  in  using  global  (bi-global)  modes  for  the  analysis  of  a  variety  of  complex 
(and  “simple” )  flows  led  to  their  applications  in  flow  control.  The  progress  in  computational 
capabilities  brings  this  advanced  technique  to  common  practice  in  studies  of  flow  pertur¬ 
bations  [The03,  EGOS,  AEGHOS,  RT08,  BABHOO,  BHHS09].  However,  the  formulation  of 
global  eigenvalue  problems  is  accompanied  by  some  uncertainties  in  the  choice  of  boundary 
conditions  over  a  2D  domain.  The  choice  of  boundary  conditions  has  a  heuristic  nature,  and 
this  provokes  questions  regarding  the  suitable  formulation  of  the  eigenvalue  problems.  In 
this  note,  I  would  like  to  suggest  a  simple  model  that  can  help  us  understand  the  effect  of 
the  upstream  and  downstream  boundary  conditions  on  the  eigenvalues  and  eigenfunctions, 
and  possible  limitations  of  the  global  formulation.  Originally,  my  idea  was  to  link  the  global 
stability  analysis  of  simple  flows  (such  as  a  flat  plate  boundary  layer)  with  the  conventional 
quasi-parallel  analysis  based  on  the  Orr-Sommerfeld  equation.  The  analysis  of  a  weakly  non¬ 
parallel  flow  has  been  carried  out  in  [AEGHOS].  A  comprehensive  analysis  of  the  upstream 
and  downstream  boundary  conditions’  effect  on  the  spectrum  of  parallel  boundary-layer  flow 
was  carried  out  in  [RT08].  These  explorations  shed  some  light  onto  the  representation  of 
discrete  TS  modes  within  the  “box  formulation.”  However,  because  the  underlying  complex¬ 
ity  of  the  2D  eigenvalue  problem  demands  signihcant  computational  power,  there  is  a  lack 
of  understanding  of  the  qualitative  features  of  this  advanced  formulation.  As  we  know,  the 
initial-value  problem  in  a  parallel  boundary  layer  flow  can  be  presented  as  an  expansion  into 
normal  modes  of  discrete  and  continuous  spectra  [Gus79,  SG81].  What  can  happen  when 
we  approach  the  problem  using  the  “box  formulation”  ?  What  can  we  get  rid  of  completely, 
and  what  are  possible  constraints  on  the  perturbations  under  consideration? 
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6.2  A  model 


To  make  the  model  as  transparent  as  possible,  I  suggest  beginning  with  the  initial-value 
problem  in  a  uniform  incompressible  flow  (f/oo  =  const).  The  solution  of  this  problem  can 
be  found  analytically  using  the  Laplace  transform  in  time  and  Fourier  transform  in  x-  and 
y-  directions.  There  are  no  discrete  normal  modes,  and  the  solution  can  be  presented  as  an 
expansion  into  the  vorticity  modes  of  the  continuous  spectrum.  For  a  given  streamwise  wave 
number,  a,  and  transversal  wave  number,  k,  the  vorticity  modes  have  complex  frequency 
uj  =  a  —  i  {k^  +  a^)  / R,  where  R  is  the  Reynolds  number  based  on  the  free  stream  velocity, 
Uoo,  and  on  the  appropriate  length  scale  L.  Let’s  address  the  known  problem  using  the 
“box  formulation”  for  a  2D  domain:  x  =  [—Lx,Lx]  ,y  =  [—Ly,Ly].  It  is  still  a  challenging 
problem  if  we  use  advanced  computational  tools  to  shed  some  qualitative  light  onto  the  effect 
of  the  boundary  conditions  on  the  2D  normal  modes.  In  order  to  elucidate  the  features  of 
using  the  “box  formulation”,  let’s  consider  a  case  when  Ly  3>  L^.  In  the  limit  of  the  “tall” 
box,  we  arrive  at  a  consideration  of  the  strip  x  =  [—1,  l],y  =  (— cxd,  -|-cx))  is  the  length 
scale).  The  limit  allows  Fourier  transform  in  y  similarly  to  analysis  in  [Tum03],  in  which  the 
strip  of  uniform  flow  with  inhomogeneous  upstream  and  downstream  boundary  conditions 
for  perturbations  was  considered.  Here,  we  use  the  same  model  of  perturbations  in  a  strip  of 
uniform  flow,  but  we  formulate  homogeneous  boundary  conditions  at  a;  =  —  1,  and  a;  =  -|-1. 

The  linearized  Navier-Stokes  equations  can  be  written  in  the  matrix-vector  form: 

dA 

ikA  -  iuHtA  =  HiA  +  H2— , 

dx 

rp 

where  H*,  Hi ,  and  H2  are  4x4  matrices.  Vector  A  =  (m,  p,  n,  D)  is  comprised  of  the  pressure 
perturbation  p,x-  and  y-  velocity  components  u  and  n,  respectively;  D  =  du/dy  —  dv/dx. 

The  solution  of  the  system  of  ODEs  with  constant  coefficients  can  be  written  in  the  form: 

4 

A  =  ^  CjZj  exp  (A^a;),  where 
i=i 
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Ai  =  -k,  A2  =  kM  =  2  ,  As  =  -  (  R  +  V4FTR^^^4L^)  , 
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The  hrst  and  the  second  fundamental  solutions  can  be  interpreted  as  “pressure  waves”  (they 
have  zero  vorticity  and  non-zero  pressure).  The  third  and  the  fourth  fundamental  solutions 
can  be  interpreted  as  “vorticity  waves”  (they  have  zero  pressure  and  non-zero  vorticity).  The 
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coefficients  Cj  are  determined  from  the  boundary  conditions.  The  box  formulation  delivers 
the  pressure  perturbation.  The  downstream  and  upstream  pressure  waves  are  localized  within 
the  layers  having  an  order  of  magnitude  ~  1/  |A;|.  In  order  to  have  a  reasonable  approximation 
of  the  unbounded  domain,  one  has  to  consider  perturbations  having  the  characteristic  scale 
/y  1,  and  to  keep  in  mind  that  the  solution  is  affected  by  the  boundary  conditions  within 
the  layers  x  ~  [—1,  —1  +  ly]  and  a;  ~  [1  —  ly,  1].  At  hnite  k  and  1,  one  can  hnd: 


X3  ^  ioo  — 


R 


A4  ~  i?  —  ioo. 


The  third  fundamental  solution  represents  a  vorticity  mode  propagating  with  the  flow  and 
slowly  decaying  due  to  the  viscosity.  The  forth  fundamental  solution  represents  a  vorticity 
perturbation  stemming  from  the  right  boundary  (a;  =  1)  and  localized  within  the  layer  of 
thickness  O  {1/R). 


6.3  An  example  of  “global”  modes  with  the  Dirichlet 
boundary  conditions 

Consider  an  example  of  the  boundary  conditions: 

u  (—1)  =  V  (—1)  =0;  u{l)  =  V  (1)  =  0. 

Homogeneous  boundary  conditions  lead  to  a  dispersion  relation  ui  =  u  {R,  k).  The  dispersion 
relation  can  be  written  in  the  form: 

A  (A:,  R,  u)  =  -2A;V4A;2  +  R^  -  4ii?cuCosh  {R)  + 

A;V4A;2  +  R^  -  4ii?o;Cosh  {2k  -  V4k^  +  R^  -  4iRu?j  + 

kV4k^  +  _  AiRuCosh  (2k  +  V4A;2  +  R‘^  -  AiRu^j  + 

2  (— 2/c^  +  iRjj)  Sinh  {2k)  Sinh  ^  V Ak‘^  +  R^  —  AiRu^  =  0 

The  eigenvalues  can  be  found  as  roots  of  A  (/c,  i?,  ui)  =  0.  One  can  also  approach  the  problem 
directly  by  using  the  spectral  collocation  method  (SCM)  for  the  system  of  ODEs.  Figure 
6.1  illustrates  the  spectrum  when  parameter  k  is  varied  (obtained  using  Mathematica)  and 
results  obtained  using  SCM  at  k  =  2  and  5. 

It  is  worth  mentioning  that  the  lines  in  Fig.  6.1  are  skewed  at  their  tips  similarly  to 
global  spectra  [see  Fig.  2  in  [AEGHOS]].  In  the  unbounded  domain,  the  continuous  spectrum 
represented  by  the  vorticity  modes  corresponds  to  a  straight  line  in  the  complex  plane  u.  In 
the  case  of  the  box  formulation,  the  tips  of  the  modes  correspond  to  small  k,  and  they  are 

affected  by  the  boundary  conditions.  Figure  6.2  illustrates  streamwise  velocity  perturbation 

in  mode  2,  k  =  5.  It  is  tempting  to  say  that  the  temporal  global  modes  also  reflect  the 
spatial  amplihcation,  but  such  a  statement  apparently  is  irrelevant  to  the  physical  problem, 
and  is  a  pure  effect  of  the  downstream  boundary  condition. 
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Figure  6.1:  Comparison  of  SCM 
results  with  results  obtained  using 
Mathematica. 


Figure  6.2:  Streamwise  velocity  per¬ 
turbation  in  mode  2,  k  =  5. 


6.4  Preliminary  conclusions 

1.  Assuming  that  the  2D  eigenvalue  problem  leads  to  a  complete  set  of  eigenfunctions, 
we  should  not  discuss  how  “wrong”  or  “correct”  the  boundary  conditions  are  in  the 
“box  formulation.”  The  discussion  should  address  the  appropriateness  of  the  set  of 
eigenfunctions  for  the  specihc  problem. 

2.  The  “box  formulation”  introduces  artifacts  such  as  upstream  pressure  and  vorticity 
waves.  There  is  a  boundary  layer  within  the  downstream  boundary  of  the  domain.  The 
pressure  wave  penetrates  upstream  at  a  distance  of  approximately  the  characteristic 
scale  of  the  perturbation  in  the  y  direction.  In  order  to  have  the  physical  picture 
unaffected  by  the  boundary  conditions,  size  of  the  box  should  be  chosen  depending  on 
the  typical  scales  of  the  perturbations  under  consideration. 

3.  The  spatial  growth  of  velocity  perturbation  within  a  global  mode  does  not  necessarily 
mean  that  the  mode  captures  a  spatial  amplihcation  associated  with  the  convective 
instability. 
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Chapter  7 
Conclusions 


A  comprehensive  study  of  stability  of  hypersonic  boundary  layers  has  been  carried  out  under 
the  support  from  the  Air  Force  Office  of  Scientihc  Research,  USAF  under  grant  FA9550-08- 
1-0332  monitored  by  Dr.  J.  D.  Schmisseur. 


The  main  results  of  the  project: 

•  The  multimode  decomposition  technique  may  serve  as  a  tool  for  gaining  insight  into 
the  flow  dynamics  in  the  presence  of  perturbations  belonging  to  different  modes.  In  the 
past,  one  could  compare  DNS  results  with  theoretical  prediction  for  the  unstable  mode 
only  far  downstream  from  an  actuator  where  the  unstable  mode  dominates  the  total 
signal.  Using  the  biorthogonal  eigenfunction  system,  one  can  compare  DNS  results 
with  theoretical  predictions  for  the  unstable  and  stable  modes  in  the  vicinity  of  the 
actuator  as  well. 

•  In  Ref.  [TWZ07]  and  in  the  present  work,  we  have  found  that  the  multimode  decom¬ 
position  requires  a  more  elaborate  analysis  within  the  point  of  synchronism  of  mode 
F  with  the  continuous  spectrum.  Therefore,  an  extension  of  the  theoretical  model  of 
Ref.  [FKOl]  to  the  case  of  continuous  spectrum  is  required. 

•  The  structure  of  the  discrete  spectrum  in  stability  analysis  of  high-speed  boundary 
layers  depends  on  basic  flow  parameters  such  as  the  Mach  number,  Prandtl  number, 
Reynolds  number  etc.  The  spectrum  may  have  two  unstable  modes  that  are  easily 
associated  with  inviscid  instabilities  of  Macks  hrst  and  second  modes.  With  another 
choice  of  the  basic  parameters,  the  spectrum  may  have  only  one  unstable  mode  having 
two  maxima  of  the  growth  rate. 

•  Terminology  introduced  by  Mack  [Mac69]  using  modes  1  and  2  (the  well-known  hrst  and 
second  modes)  is  inconsistent  with  a  proper  mathematical  treatment  of  normal  modes. 
Macks  family  of  solutions  is  actually  what  should  be  called  a  normal  mode,  whereas 
his  dehnitions  of  hrst-mode  amplihed  solutions  and  second- mode  amplihed/damped 
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solutions  lead  to  confusion  when  normal  modes  are  used  to  represent  a  solution  of  the 
linearized  Navier-Stokes  equations. 

•  The  considered  examples  of  boundary  layer  flows  in  binary  mixtures  of  oxygen  and 
nitrogen  indicate  that  the  real  gas  effects  stabilize  the  hrst  mode  and  destabilize  the 
second  one  in  the  inviscid  limit  of  the  stability  equations  (see  [FTIO]  regarding  the 
terminology).  The  result  is  in  agreement  with  [MA91]  where  this  observation  was  made 
for  boundary  layer  flow  for  air  in  chemical  and  thermal  equilibrium  with  the  stability 
equations  including  the  effect  of  viscosity,  diffusivity,  and  thermal  conductivity.  Real 
gas  effects  lead  to  a  decrease  in  wall  temperature  and  can  have  an  impact  on  the 
perturbations  similar  to  the  wall  cooling  effect  in  boundary  layers  of  non-reacting, 
calorically  perfect  gases.  However,  such  an  interpretation  of  the  wall  temperature  effect 
might  be  misleading  when  binary  mixtures  of  nitrogen  and  oxygen  are  compared.  In 
the  considered  example,  nitrogen  does  not  experience  signihcant  dissociation,  whereas 
oxygen  goes  through  dissociation.  The  difference  leads  to  qualitative  differences  in  the 
specihc  heats,  and  the  second  mode  in  nitrogen  is  more  unstable  than  the  second  mode 
in  oxygen  in  spite  of  the  higher  wall  temperature. 

•  The  “box  formulation”  used  in  global  (bi-global)  stability  analysis  introduces  artifacts 
such  as  upstream  pressure  and  vorticity  waves.  There  is  a  boundary  layer  within 
the  downstream  boundary  of  the  domain.  The  pressure  wave  penetrates  upstream  at  a 
distance  of  approximately  the  characteristic  scale  of  the  perturbation  in  the  y  direction. 
In  order  to  have  the  physical  picture  unaffected  by  the  boundary  conditions,  size  of 
the  box  should  be  chosen  depending  on  the  typical  scales  of  the  perturbations  under 
consideration. 
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Appendix  A 

Appendices  to  Chapter  2 


A.l  Comparison  of  DNS  mean  velocity  and  tempera¬ 
ture  profiles  with  the  self-similar  solution.  Com¬ 
parison  of  eigenvalues  o  =  +  iai  obtained  using 

DNS  and  self-similar  profiles. 
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Figure  A.l:  Self-similar  solution  (solid  line)  and  DNS  results  (symbols)  for  the  mean  velocity 
profile  U  {y) . 
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Figure  A. 2:  Derivative  dU/dy.  solid  line  -  self-similar  solution,  symbols  -  DNS. 


y/L 

Figure  A. 3:  Derivative  d^U/dy^\  solid  line  -  self-similar  solution,  symbols  -  DNS. 
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Figure  A. 4:  Mean  temperature  profile  T{y)\  solid  line  -  self-similar  solntion,  symbols  -  DNS. 


Figure  A. 5:  Derivative  dT/dy.  solid  line  -  self-similar  solution,  symbols  -  DNS. 
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Figure  A. 6:  Derivative  (PT/dy"^-.  solid  line  -  self-similar  solution,  symbols  -  DNS. 


Figure  A. 7:  Comparison  of  ar  obtained  using  DNS  and  self-similar  velocity  and  temperature 
profiles.  Solid  line  -  self-similar  profiles;  symbols  -  DNS  profiles. 
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Figure  A. 8:  Comparison  of  ai  obtained  using  DNS  and  self-similar  velocity  and  temperature 
profiles.  Solid  line  -  self-similar  profiles;  symbols  -  DNS  profiles. 
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